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I.  INTRODUCTION 


The  principal  kill  mechanism  of  target  structures  in  an  exoatmos- 
pheric  nuclear  environment  is  most  likely  to  be  x-ray  deposition  and  the 
accompanying  material  response  due  to  heating.  The  various  material 
responses  can  consist  of  melting,  vaporization,  shock  wave  propagation, 
and  spallation. 

The  physical  phenomenology  of  the  material  response  involved  can  be 
related  as  follows.  Initially  the  external  surface  (front  face),  assuming 
slab  geometry,  is  heated  almost  instantaneously  to  high  temperatures 
causing  material  to  melt  and  blow-off  in  the  form  of  vapor,  liquid-vapor, 
liquid,  or  solid-liquid.  The  impulsive  load  due  to  the  blow-off  creates 
a stress  wave  that  propagates  into  the  interior  of  the  target  toward  the 
bottom  surface  (back  face) . Simultaneously  the  material  below  the  top 
surface  is  heated,  inducing  a compressive  stress  which  propagates  to  the 
bottom  surface  (back  face)  and  reflects  as  a rarefaction  tensile  stress 
wave.  At  some  point  near  the  bottom  surface  this  rarefaction  tensile 
stress  wave  overpowers  the  compressive  stress  wave  moving  toward  the 
bottom  surface,  and  the  high  tension  produced  may  cause  spalling. 

The  consequence  of  the  material  response  is  that  it  leaves  the 
target  in  a weakened  condition.  In  order  to  have  a capability  to 

predict  damage  to  the  structure,  RIP^-*  ^ was  coupled  to  PETROS  3^. 

The  coupling  was  done  on  the  premise  that  calculations  of  material 
response  (which  normally  cover  a period  of  a few  microseconds,  at  most) 
and  predictions  of  structural  response  (which  generally  involve  response 
durations  on  the  order  of  milliseconds)  can  be  performed  consecutively 
rather  than  requiring  a simultaneous  treatment. 


H.  Fisher,  G.  A.  Lane,  and  R.  A.  Cecil,  "RIF,  A One-Dimensional 
Material  Response  Code  - User's  Guide,"  Systems,  Science  and  Software, 
Report  3SR-751-I,  September  1972. 

^R.  H.  Fisher,  G.  A.  Lane,  and  R.  A.  Cecil,  "RIP,  A One-Dimensional 
Material  Response  Code  - Code  Reference  Manual, " Systems,  Science  and 
Software^  Report  2SR-7 Sl-II ^ September  1972, 

R.  A.  Kruger^  ^^The  Material  Response  Program^  RIP^  " Systems^  Scvence 
and  Softu)are^  Report  3SR-120^  December  1969, 

^RlPy  A One-Dimensional  Material  Response  Code,  Systems ^ Science  and 
Softuare^  Contract  No,  SSS-R-7 2-1324, 

^S,  Atlurn^^  E,  A,  Witmer^  J,  W,  Leechj,  and  L,  Morino^  '^PETROS  3:  A Finite 

Difference  Method  and  Program  for  the  Calculation  of  Large  Elastic-Plastic 
Dynamically  Induced  Deformations  of  Multilayer  Variable  - Thickness 
Shells y " U,S,  Army  Ballistic  Research  Laboratories ^ Contract  Report 
No,  60 y November  1971  (AD  H-890200L) , 
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This  report  presents  the  calculational  procedure  used  to  predict 
the  structural  response  of  a 6061-T6  aluminum  flat  plate.  The  input 
and  results  of  RIP  calculations  are  presented  along  with  a comparison 
of  two  structural  response  calculations  utilizing  PETROS  3. 

Appendix  A contains  the  code  modifications  required  in  the  RIP  code 
to  permit  the  appropriate  data  to  be  extracted  for  use  as  input  to  the 
PETROS  3 structural  response  calculations.  In  order  to  expedite  analysis 
of  structural  response  output  data,  a plotting  program  was  developed. 

A description  of  this  program  is  presented  in  Appendix  B along  with  a 
listing.  The  BRL  version  of  PETROS  3 required  some  minimal  modifications 
and  these  are  described  in  Appendix  C.  Appendix  D contains  temperature 
variation  data  related  to  temperature  changes  which  occurred  during  the 
structural  response  calculations. 


I I . CALCULATIONAL  PROCEDURE 

A physical  description  of  the  problem  is  portrayed  diagrammatical ly 
in  Figure  1.  The  target  is  a flat  rectangular  plate  which  is  partitioned 
into  numerous  rectangular  columns  of  material  made  up  of  smaller  volumes 
of  material  called  mass  zones.  The  purpose  of  this  geometry  is  the 
accommodation  of  the  fact  that  the  RIP  code  is  one-dimensional.  Thus  for 
the  RIP  calculations , it  is  sufficient  to  concentrate  on  a specific 
column.  As  Figure  1 indicates,  the  x-rays  impinge  on  the  top  mass  zone 
of  the  column  thus  depositing  energy  and  consequently  heating  the  material 
to  high  temperatures.  A stress  wave  is  generated  and  the  situation  shown 
in  Figure  1 develops.  Some  of  the  material  is  removed  from  the  front 
face  as  liquid,  liquid-solid,  vapor  and  vapor-liquid  (for  this  work, 
only  a phase  change  to  liquid  and  liquid-solid  takes  place).  Stress 
waves  in  the  solid  material  caused  both  by  this  material  removal  and  by 
non-uniform  heating,  after  reflection  from  the  bottom  surface  (free- 
boundary), resulted  in  sufficiently  large  tensile  stress  to  cause  spalla- 
tion . 


The  RIP  code  employs  a Lagrangian  finite  difference  method  and  com- 
putes the  material  response  due  to  rapid  heating  by  solving  simultaneously 
the  equations  of  conservation  of  mass,  momentum,  and  energy. 

The  Gray*  equation  of  state  was  used  in  conjunction  with  the  simple 
von  Mises  elastic-plastic  model,  in  which  the  yield  behavior  is  calculated 
from  the  shear  modulus  and  yield  strength.  An  energy  dependent  peak 
tensile  strength  spall  model  was  used  where  the  spall  strength  is  a 
function  of  the  melt  energy.  If  the  energy  is  greater  than  the  melt 


The  Gray  equation  of  state  was  chosen  because  temperature  is  computed 
directly . If  any  other  equation  of  state  options  provided  by  RIP  were 
used^  additional  calculations  would  have  to  be  added. 
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Figure  1.  Representation  of  the  Coupling  between  RIP  and  PETROS  3 
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energy,  the  material  is  assumed  to  have  zero  strength  and  spallation 
will  be  initiated.  If  the  energy  is  less  than  the  melt  energy,  the 
tensile  stress  level  required  to  initiate  spallation  is  assumed  to  be 
the  value  of  the  spall  strength. 

The  coupling  was  performed  by  running  the  RIP  program  as  shown  in 
Figure  2,  from  t = 0,  with  an  x-ray  radiation  source  of  10  kev  blackbody 

lx  O ^ 

spectrum,  a fluence  level  of  80  cal/cm  and  a shine  time  of  75  nanoseconds, 
to  a termination  time,  t^^p.  This  time  is  determined  by  a damping  model 

described  in  Reference  6,  page  46.  This  damping  model  requires  that  each 
piece  of  solid  material  be  brought  to  a uniform  velocity,  rather  than  to 
a static  equilibrium.  It  was  incorporated  in  RIP  so  that  after 


RIP 


PETROS  3 


0 


0 

J 


PF 


RF 


Figure  2.  Time  History  of  the  Coupling 


considerable  time  transverse  (through- thickness)  wave  motion  would  vanish 
throughout  the  solid  material  leaving  the  residual  material  and  each 
spall  fragment  with  individual  uniform  velocity.  The  damping  can  be 
activated  at  any  specified  time  step  after  the  plastic  work  and  spallation 
are  essentially  complete.  This  determination  was  arrived  at  by  analyzing 
the  results  of  a run  without  damping  as  explained  in  Section  III  of  this 
report . 

*A  fluenae  level  of  80  oal/cm^  was  selected  in  order  to  guarantee  blowoff 
on  the  front  face  and  spallation  on  the  backface. 

^J.  M.  Santiago,  H.  L.  Wisniewski,  N.  J.  Buffington,  Jr.,  "A  User’s  ^ 
Manual  for  the  REPSIL  Code,"  U.S.  Army  Ballistic  Research  Laboratorzes , 
Report  No.  1744,  October  1974.  (AD  #A003176) 
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At  RIP  termination  time,  t^^p,  the  largest  residual  solid  and  its 

physical  data,  thickness,  mean  transverse  velocity,  through-thickness 
distribution  of  temperature  and  tangential  stresses  are  transferred  to  a 
rectangular  flat  plate  to  serve  as  initial  conditions  to  the  BRL  version 
of  PETROS  3 calculations  as  shown  in  Figure  1. 

PETROS  3 also  uses  a finite  difference  method.  It  solves  the  partial 
differential  equations  of  motion  that  govern  large  deflections  and  elastic- 
plastic  (von  Mises  model,  consistent  with  RIP  calculations)  behavior  of 
thin  Kirchhoff  shells  exposed  to  impulsive/pressure  loadings. 

In  working  with  the  RIP  program,  it  was  found  that  the  equation  of 
state,  (Gray  EOS) , performs  the  calculation  as  if  the  material  were 
rapidly  heated,  while  the  elastic-plastic  (von  Mises  model)  calculation 
is  carried  out  as  if  the  material  is  cold  (room  temperature)  or  at  least 
values  of  elastic  modulus  and  yield  stress  remain  constant.  The  basis 
for  this  assumption  is  that  for  the  short  time  period  covered  by  the  RIP 
calculation,  there  is  insufficient  time  for  thermal  degradation  of  material 
properties.  Proof  of  this  remains  to  be  established  for  material  which 
has  been  heated  nearly  to  its  melting  point.  However,  due  to  this  incon- 
sistency in  RIP  between  the  equation  of  state  and  the  elastic-plastic 
calculation  it  was  felt  worthwhile  to  attempt  to  bracket  the  true  response 
by  performing  two  PETROS  3 calculations  using  (a)  undegraded  material 
parameters  (results  labeled  cold  material  properties)  and  (b)  material 
properties  associated  with  long  soak  times  at  elevated  temperatures 
(results  labeled  hot  material  properties) . 

At  time  t = 0,  the  PETROS  3 program  reads  in  the  geometry  of  the 
p 

plate  and  the  physical  data  of  the  material  from  the  RIP  solution.  Mean 
transverse  velocity  is  uniform  for  every  finite  difference  mesh  point  in 

the  Y^,  plane  (see  Figure  3).  In  the  Y^  direction,  the  thickness, 
which  is  uiiiform  for  the  plate,  is  divided  into  a number  of  discrete 
layers  at  each  mesh  point. 

Values  of  tangential  stress  and  temperature*  are  constant  for  each 

layer  in  the  Y^,  Y^  plane,  but  vary  for  each  layer.  The  number  and 
thickness  of  the  layers  depend  on  the  number  of  points  one  wishes  to 
u.se  with  the  Gaussian  integration  in  determining  force  and  moment 
resultants . 


The  temperature  distribution  throughout  the  structural  response  solution 
remains  constant,  A discussion  of  this  point  is  presented  in  Appendix  C. 
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Figure  3.  Flat  Plate  Geometry  (6061-T6  Aluminum) 
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Ill,  RIP  INPUT  AND  RESULTS 


The  input  cards  listed  in  Table  I are  the  input  parameters  required 
to  perform  the  RIP  calculation.  They  are  defined  in  Reference  1,  pages 
3 through  33.  The  parameters  for  the  Gray  equation  of  state  subroutine 
are  listed  in  Table  II. 

Figure  4 shows  the  energy  per  gm  (dose)  deposited  through  the  initial 
thickness  (5  mm)  which  was  computed  by  the  energy  deposition  and  transport 
7 8 

program  FSCATT  ■* *  . This  program  serves  as  the  x-ray  radiation  source 

for  RIP.  Figure  5 shows  a representation  of  the  column  element  of 
material  at  the  end  of  the  RIP  calculation.  It  has  gone  through  various 
phase  changes  due  to  x-ray  energy  deposition.  The  loss  of  material  due 
to  blowoff  on  the  front-face  is  in  the  form  of  liquid  and  liquid-solid 
material,  while  the  loss  of  material  on  the  back-face  is  caused  by 
spallation  in  the  solid  region. 

The  RIP  code  calculates  the  hydrostatic  pressure  P and  the  deviatoric 
stress,  SD,  (which  is  actually  the  component  of  deviatoric  stress,  , 

associated  with  the  through-thickness  direction) . The  transverse  (through- 
thickness) stress,  (see  Figure  5)  is  given*  by  = a^'  + P.  In 

order  to  determine  residual  stresses,  the  RIP  code  has  been  modified  to 
provide  for  damping  out  the  transverse  stress  waves.  Figure  6 shows  the 
distribution  of  transverse  stresses  at  t^^p.  It  may  be  seen  that  these 

stresses  are  negligible  throughout  the  thickness  of  the  solid  material. 
However,  the  corresponding  hydrostatic  pressures  shown  in  Figure  7 are 
non-negligible . 

Since  we  know  that  the  sum  of  the  deviatoric  stress  components  must 
vanish  and  that  the  tangential  stresses,  a^,  o^,  must  be  equal  for  the  one- 
dimensional RIP  analysis,  it  follows  that  = 02  " ^ ^ plotted 

in  Figure  8.  It  should  be  noted  that,  for  the  stresses  fully  damped 
out,  = 1-5  P.  These  tangential  stresses  are  the  combined  result 

of  the  inelastic  deformation  produced  by  the  transverse  stress  waves 

and  the  thermal  stresses  corresponding  to  the  temperature  distribution 

at  t„„  (see  Figure  9) . 

Kr 


^i?.  H.  Fisher  and  J.  W.  Wiehe,  "A  User's  Guide  to  the  FSCATT  Code," 
Systems,  Science  and  Software,  Report  3SR-328,  November  1970. 

^R.  H.  Fisher  and  R.  A.  Kruger,  "A  Numerical  Treatment  of  Scattering 
and  Fluorescence  in  Plane  Geometry, " Systems,  Science  and  Software, 
Report  3SR-119,  September  1969. 

*For  the  RIP  analysis,  compressive  stresses  are  considered  positive. 
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Table  I.  Input  for  RIP  Calculation 


( AMASS ) 

Bl=1.352 

( MATBND) 

B31=10C 

( MATLID) 

B41=1C3 

( NMTRLS) 

861  = 1 

( MAXCYCL  ) 

B84=995 

( MoDE ) 

B85  = 0 

( TMAX ) 

B87=6.0E-04 

( IQUNIT) 

B245=2^-1,3*0,4 

(NAPS ) 

8291=1,17 

( NORD  ) 

8296=10,11 

( NCUMP  ) 

8313=1000 

(NPRINT) 

8315=100 

( TPR  INT) 

8316=1. 92396E-06 

( niSCPT) 

B1329=12A 

X-RAY  LOADING  ON  6061-T6  ALUMINUM 
( ICCMBZ ) X1051  = l 

(SPALL  CRITERIA)  X1052=I 

X1001=3 
B182=1C.0 
B185=-80.0 
B219=0. 

B221=.75E-07 
B227=l 

V66  1 = 2. 704,0. ,0. ,0. ,C. , 5. 2 1 G5 , 2. 75 1 62 E 1 1 , 2. 69948 E9 , 
-1.3E10, 1.2E11,3.C5E10, 1. 1F10,7. 15E9,0. ,0., 

1.338, 2. 18, 1. 68,. 666667, 8. 7E-9, 1340. ,0. ,0. ,26.98 , 
.452,.  190, 47.,!., 11. 46966-5, 3. 

(END  OF  DATA)  END 

1 4Cj  X 

^ I See  Appendix  A for  detailed  explanation. 

l.eCE-06  10.0* 

BRCCYX  -1 

0 


(NSPALL) 

( BBTEMP) 

( ENERGYl  ) 
( SSTARTl  ) 
( SSTOPl) 
(NFIT) 
(GtCSAL) 


1 

14 

.006 

1 

12 

.01 

1 

29 

. 0025 

1 

24 

.002 

1 

13 

.9795 
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Table  II.  Gray  Equation  of  State  Parameters  for  6061-T6  Aluminum 
Numerical  Value  Property 


2.704 

gm/cm^ 

9 

Ambient  density 

5.21  X 10^ 

cm/ sec 

g 

Bulk  sound  speed 

2.751  X 10^^ 

dynes/cm^ 

g 

Shear  Modulus 

2.699  X 10^ 

dynes/cm^ 

g 

Yield  strength 

-1.3  X 10^° 

dynes/ cm^ 

2 

Spall  limit 

1.2  X 10^^ 

ergs/cm 

2 

Material  has  vaporized 

3.05  X 10^° 

ergs/cm 

2 

Material  commences  to  vaporize 

1.10  X 10^° 

ergs/cm 

2 

Material  has  completed  melting 

7.15  X 10^ 

ergs/cm 

2 

Material  commences  to  melt 

1.338 

12 

Hugoniot  parameter 

2.18 

dimensionless 

1 2 

Lattice  gamma 

1.68 

^12 

a = Yo  - -5 

.6667 

Electronic  gamma^"^ 

8.7  X 10”^ 

3 2 

Mbar  cm  /mole  deg 

10 

Electronic  energy  coefficient 

1340. 

12 

Melting  temperature  parameter 

0. 

^ 2 
Mbar  cm  /gm 

12 

Energy  at  reference  state 

26.98 

gm/mole 

Atomic  weight^^ 

.452 

cm^/gm 

1 2 

Volume  at  which  EOS  are  joined 

. 190 

cm^/gm 

12 

Excluded  volume  for  vapor  phase 

47. 

3 

Mbar  (cm  /mole) 

Coefficient  of  attractive 

12 

potential  for  vapor 

1.0 

dimensionless 

T • ^ 

Join  parameter 

11.469  X lo"'^ 

3 

Mbar  cm  /mole  deg 

12 

Entropy  of  melting 

9 

B,  J.  Kohn^  ^^Compilation  of  Eugoniot  Equations  of  State^  ''  Air  Force  Weapons 
Laboratory^  April  1969  (page  19). 

^^Sensitivity  of  Material  Response  Calculations  to  the  Equation  of  State  Models 
Systems^  Science  and  Softuare^  Contract  Report  Eo.  130  (SSS-‘R~73~1910 ) ^ 
December  1973  (page  24). 

^^E.  B.  Royce^  *^Gray^  A Three-Phase  Equation  of  State  for  Metals Lawrence 
Livermore  Laboratory^  Report  UCRL-61121^  September  1971  (page  38). 
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Dynamic  Response  of  Aluminum  to  Pulsed  Energy  Deposition  in  the  Melt-Dominated 
Regimey^'  U.S.  Harry  Diamond  Laboratories y HDL-TR-1624y  July  1973  (page  11). 
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Figure  4.  10  Kev  Blackbody  Dose  Profile 
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Figure  5.  Material  Break-Up 
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Figure  6.  Transverse  Stress  Distribution 
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Figure  7,  Pressure  Distribution 
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Figure  8.  Residual  Tangential  Stress  Distribution 
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Figure  9.  Temperature  Distribution 
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DISTANCE  FROM  FRONT  FACE, 


As  previously  mentioned,  damping  is  introduced  into  the  RIP  code 
calculations  after  plastic  flow  and  spallation  have  subsided.  This  is 
illustrated  in  Figure  10,  which  presents  the  kinetic  energy  of  the  solid 
material  as  a function  of  time.  Where  there  is  more  than  one  curve  for 
a period  of  time,  the  kinetic  energies  of  the  residual  material  and  of 
one  or  more  spall  fragments  are  displayed  individually.  The  lowest 
curve  gives  the  kinetic  energy  of  the  residual  plate.  Spallation  occurs 
in  the  manner  illustrated  in  Figure  5 with  separation  and  coalescence  as 
indicated  in  Figure  10.  This  case  was  initially  run  without  damping  to 
2,400  nanoseconds  (plotted  as  dashed  lines).  By  examination  of  this 
solution  it  was  determined  that  plastic  work  and  spallation  is  completed 
by  1,800  nanoseconds.  Therefore  it  would  be  satisfactory  to  introduce 
damping  at  1,800  nanoseconds.  This  was  done  in  the  second  run  which 
self-terminated  at  4,298  nanoseconds  when  a criterion  for  smallness  of 
residual  oscillations  was  satisfied.  The  corresponding  variation  of  the 
velocity  of  the  mass  center  of  a portion  (to  be  defined)  of  the  solid 
material  is  shown  in  Figure  11.  Since  the  thickness  of  residual 
material  was  not  known  until  the  completion  of  the  calculation,  it  was 
not  practical  to  plot  the  mass  center  velocity  of  that  material  alone. 
Consequently,  Figure  11  includes  all  the  variable  mass  of  solid  material 
which  has  not  spalled  off  at  the  particular  time.  After  2,500  nanoseconds 
the  plot  represents  the  velocity  of  the  residual  solid  material. 


IV.  PETROS  3 INPUT  AND  RESULTS 

PETROS  3 code  modifications  designed  to  accept  the  RIP  output  data 
are  explained  in  Appendix  C.  The  normal  PETROS  3 input  is  described  in 
Reference  5,  pages  160  through  171.  To  take  advantage  of  symmetry  only 
a quarter  of  the  plate,  as  shown  in  Figure  3,  is  treated.  The  problem 
is  solved  by  utilizing  a 12  x 16  square  mesh  with  one  layer  of  material 
and  a thickness  of  4.093  mm  (0.16116  inches)  (see  Figure  12).  An 
impulsive  velocity  of  0.83439  m/sec  (32.85  in/sec)  was  distributed 
uniformly  over  the  surface*  of  the  plate. 

In  the  structural  shell  response  analysis,  stresses  are  calculated 
ad  discrete  layers  in  the  through-thickness  direction.  These  calculations 
can  be  performed  using  "room  temperature"  material  properties  or  thermally 

IZ 

degraded  properties.  For  the  latter  case  yield  stress  and  elastic  moduli 


Along  the  vow  of  mesh  points  adjacent  to  a clamped  boundary  the  PETROS  Z 
code  will  modify  this  initial  velocity^  for  compatibility  with  the  finite 
difference  clamped  edge  boundary  condition. 

^^Engineering  Data  for  Aluminum  Structures^  The  Aluminum  Association, 

New  York,  N.  Y.,  August  1969. 

Lipkin,  J.  C.  Swearengen,  C.  H.  Karnes,  ''Mechanical  Properties  of 
6061-T6  Aluminum  after  Very  Rapid  Heating,”  Sandia  Laboratories 
Report  SC-RR-7 2-0020,  March  1972. 
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Figure  10.  Kinetic  Energy  of  Solid  Material 
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Figure  11.  Velocity  of  Mass  Center  of  Unspalled  Solid  Material 
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Figure  12.  Finite  Difference  Mesh  for  Undeformed  Plate 
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as  function  of  temperature,  as  represented  by  Figures  13  and  14,  were 
used.  Poisson *s  ratio  was  calculated  as  a function  of  temperature  from 
the  data  of  Figure  14  by  the  known  relationship  among  the  elastic 
constants  of  an  isotropic  material: 


In  PETROS  3 a Gaussian  quadrature  integration  as  explained  in  Reference 
15,  page  41,  is  used  to  evaluate  the  inplane  stress  and  moment  resultants 
at  each  mesh  point.  For  the  results  which  follow,  values  of  stress  at 
six  Gauss  points  were  used  for  the  initial  tangential  stress  through  the 
thickness  of  the  plate.  The  values  of  both  cold  and  heated  material 
properties  used  in  these  calculations  are  listed  in  Table  III  for  each 
Gauss  point  through  the  thickness. 

The  deflection  patterns  predicted  for  cold  and  heated  material 
properties  are  compared  at  two  successive  times  in  Figures  15  and  16. 

On  the  left  of  each  isometric  plot  the  deformed  and  undeformed  cross- 
section  and  profile  of  the  plate  can  be  compared.  The  deflection  of 
the  midpoint  of  the  plate  as  a function  of  time  is  shown  in  Figure  17. 

For  both  solutions,  the  plate  initially  moves  downward  due  to  the  blow-off 
impulse  but  the  locked-in  stresses  (from  the  RIP  solution)  prevail  and 
the  plate  eventually  assumes  a buckled  upward  configuration.  Not  sur- 
prisingly, the  excursions  in  both  directions  are  greater  for  the  thermally 
degraded  properties  solution. 

1 2 

Figure  18  shows  elongational  strain  versus  time  in  the  Y and  Y 
directions  (see  Figure  12)  at  the  center  of  the  plate  on  the  upper  and 
lower  surfaces  for  cold  and  heated  material  properties.  In  comparing 
these  strains,  one  sees  that  the  largest  strain  component  occurs  on  the 

lower  surface  in  Y^  direction. 

Figure  19  shows  the  energy  balance  information  produced  by  the 
PETROS  3 code  for  the  cases  under  consideration.  These  results  are 
useful  both  for  detection  of  possible  numerical  instabilities  and  for 
determining  when  the  solution  may  be  terminated.  The  information  for 
the  cold  material  properties  shows  that  the  response  was  principally 
elastic  with  very  little  plastic  work  having  taken  place  during  the 
structural  response  phase.  By  contrast,  a large  amount  of  plastic  work 
is  evident  for  the  solution  using  thermally  degraded  material  properties.  t- 

This  plastic  work  can  be  interpreted  as  being  composed  of  two  parts. 


L.  Morino^  J.  W,  Leech  and  E,  A.  Witmer^  ^'PETROS  2:  A New  Finite- 

Difference  Method  and  Program  for  the  Calculation  of  Large  Elastic- 
Plastic  Dynaxrrically- Induced  Deformations  of  General  Thin  Shells^  " 

BEL  CP  12  (MIT-ASRL  TR  152-1)^  December  1969.  (In  two  parts:  AD  708773 
and  AD  708774). 
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Figure  13.  Yield  Stress  Dependence  on  Temperature 
for  6061-T6  Aluminum  Alloy 
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From  front  face. 
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CYCLE  600 

TIME  = 2760  microseconds 
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Figure  15.  Deformed  Reference  Surface  at  Time  Step  600 
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CYCLE  1200 

TIME  = 5520  microseconds 
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Figure  16.  Deformed  Reference  Surface  at  Time  Step  1200 
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Figure  17.  Transverse  Displacement  of  Midpoint  of  Plate 
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Figure  18.  Surface  Strain  on  the  Top  and  Bottom  Surface 
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interpreted  as  being  composed  of  two  parts,  separated  by  the  dashed  line 
in  Figure  19(b).  Recalling  that  the  initial  stresses  for  the  PETROS  3 
analysis  were  obtained  from  the  RIP  solution  at  t^^p  (which  used  cold 

material  properties  for  stress  evaluation) , it  is  easy  to  see  that  the 
PETROS  3 analysis  would  entail  large  changes  in  stress  during  the  first 
time  step  when  thermally  degraded  properties  are  used.  These  essentially 
discontinuous  stress  changes  entail  a significant  reduction  in  strain 
energy  during  the  first  cycle  of  PETROS  3 calculations  and  account  for 
the  area  in  Figure  19(b)  labeled  ’^mismatch’*  plastic  work.  This  non- 
physical discontinuous  change  could  have  been  avoided  if  the  RIP  code 
had  provision  for  using  thermally  degraded  material  strength  properties. 
The  additional  plastic  work  associated  with  the  structural  response 
phase,  using  thermally  weakened  material  characteristics  is  essentially 
complete  by  the  termination  of  the  solution.  The  increase  in  external 
work  is  due  to  the  change  in  Poisson’s  ratio  from  cold  to  heated 
material  properties. 

It  should  be  noted  that  the  PETROS  3 calculations  were  performed 
using  the  temperature  distribution  of  Figure  9 as  though  the  temperature 
were  invariant  £or  the  duration  of  the  structural  response  prediction. 

A discussion  of  this  point  is  presented  in  Appendix  D. 


V.  CONCLUSION 

The  present  work  illustrates  a procedure  for  predicting  the 
structural  response  of  target  structures  after  intense,  rapid  heating 
due  to  x-ray  radiation  deposition.  It  is  evident  from  comparing  the 
solutions  of  the  undegraded  material  properties  and  the  degraded 
material  properties  that  the  results  are  significantly  different.  This 
indicates  that  precise  modeling  of  structural  response  in  rapid  heating 
situations  requires  an  accurate  definition  of  degraded  material  properties 
under  such  conditions. 

The  methodology  described  is  not  limited  to  cases  solved  here. 

Other  initial  geometries  with  multilayer  hard  bonded  material  could  be 
analyzed  in  a similar  manner. 
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APPENDIX  A 


CODING  CHANGES  TO  RIP 

Two  types  of  changes  to  the  RIP  program  listed  in  Reference  2 were 
needed  in  order  to  couple  RIP  and  PETROS  3.  The  first  change  was  de- 
signed to  extract  data  to  assist  in  interpreting  the  RIP  output  needed 
as  PETROS  3 input.  The  second  was  to  extract  temperature  from  the  Gray 
equations  of  state,  and  change  the  momentum  editing  flags  to  agree  with 
the  Gray  equation  of  state  editing  flags,  which  determine  within  which 
region  the  calculation  point  lies. 

A list  is  provided  with  a description  of  variables  used  in  sub- 
routine DAMP,  PDATA,  and  PETINP.  They  are  constructed  to  speed  the 
process  of  determining  the  solution  and  provide  additional  output  not 
contained  in  the  normal  RIP  output  format.  Following  this  list  is  a 
brief  explanation  and  FORTRAN  listing  of  the  subroutines. 
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CHANGES  TO  THE  RIP  PROGRAM  TO  OBTAIN  PLOTTING  DATA  FOR  THE  BRL  R 
PROGRAM 


1.  FOLLOWING  CARD  RIP1270  ADD  THE  FOLLOWING 

IFINRESRT  .EQ.  OICALL  PDATA( 1 » 

2.  FOLLOWING  CARD  RIP1310  ADD  THE  FOLLOWING 

CALL  PDATAin 

3.  FOLLOWING  CARD  RIP1370  ADD  THE  FOLLOWING  STATEMENTS 

CALL  EDIT 
CALL  TDELT 

A.  FOLLOWING  CARD  RIP2190  ADD  THE  FOLLOWING  STATEMENTS 
ICO  CALL  PDATA(2» 

CALL  EDIT 
CALL  TDELT 
CALL  PDATA(5> 


REMOVE  CARD  RIP2200  WHICH  STATES 
180  CALL  EDIT 

5.  FOLLOWING  CARD  EDI  950  ADD  THE  FOLLOWING 

CALL  PDATA(4» 

6.  FOLLOWING  CARO  EDIIAIO  ADO  THE  FOLLOWING 

CALL  PETINP 

IFINSEDIT  .EQ.  0)CALL  PDATA(3) 

CALL  EXIT 

REMOVE  CARD  EDIIA2C  WHICH  STATES 
IFINSEDir  .EC.  0)CALL  EXIT 


CHANGES  TO  THE  MOMENTUM  ED  I T I NG  FLAGS  IN  SUBROUTINE  ENCALE 

1.  CODING  WHICH  FOLLOWS  LINE 

TEMPI  1 ) = TIME-DTNP1 

ADD  L INE 

GOTO  (U2,1A2,1A6,1A6»,IE0S 

AND  CHANGE  THE  STATEMENT  WHICH  STATES 

IF(EIJ)  .GE.  ESTCONI 10,N))M20NE(M,1)=J 

TO  THE  FOLLOWING 

1A2  IFIEIJI  .GE.  ESTCCN(10,NnH20NE(M,l»=J 

2.  CODING  WHICH  FOLLOWS  LINE 

IFISIJ)  .GE.  ESTCONI 13, N) IM20NE (M,A)=J 

ACD  L INES 
GOTO  150 


FOR 

GRAY  E.O.S.  ONLY 

1F( 

IVFLAGI 

[ J) 

.EO. 

35)HZOME(P, \ )=J 

IF  ( 

IVFLAGI 

[ J) 

• E0« 

33)HZ0ME(H,2)=J 

IF( 

IVFLAGI 

[ J) 

.EQ. 

32)MZ0NE(K,3)=J 

IF( 

IVFLAGI 

[ J) 

• EQ. 

31  )HZOrjE(  K,4)=J 

P PLOTTING 


ENC1A40 


ENC1460 


ENC1490 
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CHANGES  FOR  EXTRACTING  TEMPERATURE  FROM  THE  GRAY  EQUATION  OF  STATE 
SUHROUTINE  ECS3 

1.  CODING  WHICH  FOLLOWS  LINE 

. U( 100C)fX( ICOO) »8LANK(39) 

ADD  LINE 

. COMMON  /ECSMAT/GTEMPdOOO) 

2.  CODING  WHICH  FOLLOWS  LINE 

P( J)=TEM(A)*TEHP( 5)-TEMP(7)fTEMP(9)*TEHP( 10) 

ADD  LINE 

TEMPI 5)=TEMP<8) 

3.  CODING  WHICH  FOLLOWS  LINE 

P( J)=l.E^l2*P( J) 

ACD  L IME 

GTEMPI J)=TEMP(5) 


CHANGES  FOR  PRINTING  TEMPERATURE  IN  SUBROUTINE  MECIT 

U CODING  WHICH  FOLLOWS  LINE 

. U(  ICCC)  fXUCOO)  »BLANK(39) 

ACD  LINE 

COMMON  /ECSMAT/GTEMPdOOO) 

2.  CHANGE  THE  STATEMENT  WHICH  STATES 

PRINT  13C»  Jf  SMASSfXU)  ,U(J)  »RHO(J)  ,Q(J)»PTfP(J)9E(J)tC(J) 

• IVFLAG(J),IBFLAG(J),J 
TO  THE  FOLLOWING 

PRINT  I30»J«SMASS»X(  J)  fU(  J)  »RHO(J),SD(J)  fS(J)tPU)fE(J)  yGTEMPI  J)t 

• I VFLAG(J) tlQFLAGI J) , J 

3.  CHANGE  FORMAT  ICO  TO  THE  FOLLOWING 

ICO  FORMAT! IHI, AX, 12A6,5X,IIHCYCLE  NO,  = I 5 , 2X , 6HT I ME  = IPE 1 3 • 5// AX, 
•IHJf AX.5HSMASS,8X, IHX, 10X,IHU,10X,3HRH0,7X,2HSD,9X,IHS, I IX, 

•IHP, lOX, IHE, 10X,IHT,8X,2HIV,2X,2HIB,AX, IHJ/) 

A.  CHANGE  FORMAT  lAO  TO  THE  FOLLOWING 

I AO  FORMAT! IH 1 , 3X , IHJ , AX , 5HSMASS , 8X , IHX, I OX , LHU , lOX , 3HRH0 , 7X, 2HSD,9X, 
• IHS, IIX, IHP, 10X,IHE, 10X,IHT,3X,2HIV,2X,2HI B,AX, IHJ/ ) 


ESC  210 


ESC2190 


ESC22A0 


MOT  210 


HOT  390 
MOT  AOO 
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Subroutine  PDATA,  DAMP,  and  PETINP  Local  Variables 


Name  Definition 


CMTOIN 

Constant  to  convert  centimeters  to  inches. 

CINl 

Kinetic  energy  of  the  vapor-liquid  region. 

CIN2 

Kinetic  energy  of  the  liquid  region. 

CIN3 

Kinetic  energy  of  the  liquid  solid  region. 

CINET 

Total  kinetic  energy  of  all  regions. 

CINS(20) 

Kinetic  energy  of  each  solid  piece  KF,  after  spallation. 

DCPSI 

2 

Constant  to  convert  dynes/cm  to  psi. 

DFACT 

Minimum  average  transverse  stress  used  to  control  the  ter 
mination  of  the  problem  during  damping  (Input  card  4) . 

DST0P(20) 

Control  if  DSTOP(I)  = 1 solid  piece  has  been  damped, 
DSTOP(I)  = 0 continue  damping. 

EERGY 

Internal  energy  plus  kinetic  energy  for  all  regions. 

GCTOLB 

2 2 4 

Constant  to  convert  g/cm  to  lb-sec  /in.  . 

GTEMP(J)  Temperature  of  Zone  J (K) 
J = 1,  1000 


I 

Counter. 

IFIRST 

Control  variable  governing  mode  of  operations. 

IFLAG 

Control  variable  (see  description  of  variables  in 
Appendix  B) . 

J 

Counter. 

JK 

Counter  for  indexing  array  JSS(JK) . 

JS 

Initial  count  used  in  DO  statement  in  subroutine  DAMP. 

JSN 

Final  count  used  in  DO  statement  in  subroutine  DAMP. 

JSS(40) 

Zone  numbers  corresponding  to  the  first  and  last  zone  of 
each  solid  piece. 

JSS(JK)  = First  zone  of  the  KF  solid  piece. 

JSS(JK  + 1)  = Final  zone  of  the  KF  solid  piece. 
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KEY 

Control  variable  governing  mode  o£  operations. 

KF 

Number  o£  solid  pieces  a£ter  spallation. 

KINETl (20) 

Kinetic  energy  at  time  t - At  £or  the  KF  solid  piece. 

KINET2(20) 

Kinetic  energy  at  time  t £or  each  solid  piece  KF. 

LINK 

Control  number  used  in  subroutine  PDATA. 

LMAXP 

Number  o£  times  at  which  pro£ile  plots  are  drawn  (Input 
card  2 needed  £or  plotting) . 

LTP 

Control  number. 

NCOM(20) 

Control  number  to  check  i£  two  solid  pieces  have  recombined 
a£ter  they  had  previously  spalled. 

NP. 

Number  o£  time  steps  at  which  pro£ile  plots  are  drawn 
(Input  card  1 needed  £or  plotting) . 

NPOLT 

Plotting  tape  unit  number. 

NT 

Counter . 

NTIME(50) 

Time  steps  at  which  pro£ile  plots  are  obtained  (Input  card  1) 

SMASS 

Sum  o£  the  masses  in  the  solid  region. 

SUMKIN(20) 

Sum  o£  kinetic  energies  o£  the  solid  region. 

TANGS 

Tangential  stress. 

TDAMP 

Time  at  which  damping  starts  (Input  card  4) . 

TEMF 

Temperature  in  (°F) . 

THICK 

Thickness. 

THICKI 

Spatial  coordinate  o£  the  £irst  zone  o£  the  KF  solid  piece. 

TPLOT(50) 

Times  at  which  pro£ile  plots  are  obtained  (Input  card  3) . 

TRANS 

Average  transverse  stress. 

VCEN(20) 

Array  containing  the  velocity  o£  the  mass  center  o£  each 
solid  piece  KF. 

VCENT 

Same  as  VCEN(20) . 

VCENTR 

Velocity  o£  the  mass  center  o£  the  £irst  KF. 
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Explanation  of  Subroutines 


PDATA  Along  with  the  RIP  input  data,  Table  I,  the  following  input  data 

is  needed  to  obtain  data  on  tape.  Also  this  subroutine  calculates 
and  stores  the  data  on  tape  required  by  the  BRL-RIP  plotting  pro- 
gram. Appendix  B describes  and  lists  the  program. 


Input  Data 

Card  Variables  Format 

1 NP,  (MTIME(I),  1=1,  NP)  1615 

2 LMAXP  1615 

3 (TPLOT(I),  1=1,  LMAXP)  1615 

4 TDAMP,  DFACT  2E12.6 

Description  of  Input 

Card  1 NP  Number  of  time  steps  for  which  timewise  plots  of  trans- 

verse stress,  tangential  stress,  pressure  and  tempera- 
ture are  plotted. 

NTIME(I)  Time  steps  that  transverse  stress,  tangential  stress, 
pressure  and  temperature  are  plotted. 

Card  2 LMAXP  Number  of  times  for  which  timewise  plots  of  transverse 
stress,  tangential  stress,  pressure  and  temperature  are 
plotted. 


Card  3 TPLOT(I)  Times  that  transverse  stress,  tangential  stress, 
pressure  and  temperature  are  plotted. 


Card  4 TDAMP  Time  at  which  damping  starts.  See  Reference  9,  page  59, 
for  a more  detailed  explanation  of  TDAMP. 


DFACT  Minimum  average  transverse  stress  value  used  to  control 
the  termination  of  the  problem  during  damping. 


DAMP  As  already  mentioned  in  the  Calculational  Procedure,  this  sub- 
routine is  patterned  on  subroutine  DAMP,  Reference  6,  page  46, 
except  that  no  energy  is  dissipated  by  viscous  damping,  all  the 
kinetic  energy  is  reduced  by  the  kinetic  energy  annihilation 
procedure.  This  procedure  involves  freezing  the  zone  velocity 
of  the  material  to  the  velocity  of  the  mass  center  whenever  the 
kinetic  energy  of  the  material  reaches  a maximum.  Therefore, 
the  vibratory  kinetic  energy  of  the  material  is  made  to  vanish 
and  the  total  kinetic  energy  reduced  to  the  translational  kinetic 
energy  of  the  mass  center.  Then  the  material  is  released  with 
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less  energy.  The  above  is  repeated  until  the  average  transverse 
stress  is  less  than  DFACT.  If  spallation  in  the  solid  region 
takes  place,  each  solid  portion  is  damped  separately. 

PETINP  At  the  termination  of  the  run,  this  subroutine  calculates  and 

outputs  the  physical  properties  of  the  remaining  solid  portions 
of  material  (See  Table  A-1)  which  are  needed  by  PETROS  3 as 
initial  conditions. 


47 


Table  A-1,  Rip  Output  Summarizing  the  PETROS  3 Input 


FETROS  INPUT 


SPALL  FRAGMFNT  « 1 

velocity  of  mass  center  • .32S5S6E  02  KINETIC  E • .3S6H2E  04 


transverse 

tangential 

bulk 

density 

J 

X (INCHES) 

STRESS  (PSll 

STRESS  (PSD 

MODJLUS  (PSD 

5 

.286972E-0? 

,356563E-0l 

17i926E 

05 

,99166iE 

07 

,243965E-03 

6 

,49n975E-n2 

, 177426E 

01 

181216F 

05 

, 101072E 

08 

, 246288E-03 

7 

.693n53E-02 

,2S5638F 

01 

192885E 

05 

, 101796E 

08 

, 247  166E-03 

8 

,89441 3E-02 

.251745E 

01 

17241  IE 

05 

, 1 n?03nF 

08 

,247454E-03 

9 

, 109554E-01 

.553275E 

01 

152215E 

05 

, 102100E 

08 

,24754iE-03 

10 

, 12O659E-01 

,?84710E 

01 

135297E 

05 

, 1021  lOE 

08 

, 247557E-03 

1 1 

. 1 4O764E-0 1 

,7676SnE 

01 

143511E 

05 

, 102132E 

08 

,247582E-03 

12 

, lS0866E-n 1 

.401311E 

0l 

185541E 

05 

, 102165E 

08 

,24761RE-03 

13 

, 189965E-01 

,896  145E 

0 1 

183026E 

05 

, 1 02304E 

08 

,247784E-03 

14 

,2ino5iE-oi 

,549617F. 

01 

196145E 

05 

, 1 02540E 

08 

, 248068E-03 

15 

,230li4E-Ol 

, 103269E 

02 

?068?nE 

05 

, J02753E 

08 

,248324F-03 

IS 

, 2501 5SE-01 

,5370 24 E 

01 

215400F 

05 

, 1 02945E 

08 

,246556E-03 

17 

,2701 79E-01 

• 122974E 

02 

219507E 

05 

, 1 031  14E 

08 

,248760E-03 

18 

,200186E-01 

, 439690E 

01 

218316E 

05 

, IO3261E 

08 

,246938E-03 

19 

y3ini7OE-0i 

.122163F 

02 

212850E 

05 

, 103387E 

08 

,249n92E-03 

20 

,3301 60E-01 

,591 080F 

01 

204535E 

05 

, 103498E 

08 

,249227E-03 

21 

.35O129E-01 

.813348E 

01 

19607  IF 

05 

, 103597E 

OR 

,249349E-03 

22 

, 37n08OE-01 

,9210 1 3E 

01 

18791 4E 

05 

, 103687E 

08 

,249458E-C3 

23 

,39nn40E-oi 

.355247E 

01 

180423E 

05 

, 103770E 

08 

,249560F-03 

24 

,40«983E-01 

, 1 09734F 

02 

173155E 

05 

, 1 03849E 

08 

, 249658E-03 

25 

,42991 8E-01 

, 1 19581E 

01 

186274E 

05 

, 103925E 

08 

,249750E-03 

2S 

, 44H848E-0 1 

.114176E 

02 

159608E 

05 

, 1 03995E 

08 

, 249836E-03 

27 

, 46Q766E-0 1 

.534390E 

no 

153337E 

05 

, 1 04062E 

08 

,249916E-03 

28 

,489C81E-01 

,111907E 

02 

1 46965E 

05 

, 1 04  124E 

08 

,249994E-03 

29 

,5005R9E-01 

,20n079F 

01 

14  1081F 

05 

, 1 04  183E 

08 

,250066E-03 

30 

,52049?E-0 1 

.Q51144F 

0 1 

134S17E 

05 

, 104237E 

08 

,250i33E-03 

31  . 

,54038QE-01 

,605G53E 

01 

129195E 

05 

, 1 0428OE 

08 

,250196E-03 

32 

,58928! E-01 

, 586300E 

0! 

123232E 

05 

, 1 04337E 

08 

,250255E-03 

33 

,589i6OE-0l 

.823557E 

01 

1 17596E 

OS 

, 1 04383E 

08 

,250312E-03 

34 

,60qo52E-01 

-,67723iE 

on 

1 12336E 

05 

, 1 04426E 

08 

,250365E-03 

35 

,62HQ30E-0l 

• 1 004  44F 

02 

107072E 

05 

, 1 04468E 

08 

, 25041 7E-03 

36 

,64H805E-01 

130240E 

01 

lOl 850E 

05 

,104508E 

08 

,250466E-03 

37 

,66867SE-01 

.787688E 

01 

968894E 

04 

,104547E 

08 

,2505] 3E-03 

38 

,680543E-01 

-.234712F 

01 

Q28575E 

04 

, 104584E 

08 

,250559E-03 

39 

,7084f)6E-01 

.754370F 

01 

886126E 

04 

, 1 04620E 

08 

,250603E-03 

40 

,7?8266E-01 

, 1 1 0207E 

01 

821242E 

04 

,10465!E 

08 

,25064iE-03 

41 

,74M123E-01 

,2l0b68E 

01 

6 8 5 1 2 4 K 

04 

, 104670E 

08 

,250667E-03 

42 

. 7S7978E-01 

,49?710E 

01 

515550E 

04 

, 104682E 

08 

, 250684E-03 

43 

,787831 E-0 1 

,111 253F 

01 

351117E 

04 

, 104693E 

OR 

, 250700E-03 

44 

,807683E-n 1 

, 543446F 

01 

225231E 

04 

, 1 04708E 

08 

,250720E-03 

45 

,827534E-01 

,304993E 

on 

888539E 

03 

, 10472OF 

08 

,250737E-03 

46 

,847383E-01 

,0209! 7E 

01 

643682E 

03 

, 104728E 

08 

,250750E-03 

47 

,86/23lE-ni 

,?5?6?6r-0l 

21 121  IF 

04 

, 104736F 

08 

,25076?E-03 

48 

,887079E-n  1 

,715619F 

01 

30991 4E 

04 

, 104751E 

08 

,250782E*03 

49 

,906Q24E-01 

,488688F 

01 

385713E 

04 

, 104769E 

08 

,250805F-03 

50 

,926768E-01 

,343853F 

01 

453886E 

04 

, 1 04787E 

08 

,250828E-03 

51 

,9466i0E-0l 

,61  1903E 

01 

527606E 

04 

, 1 04803E 

08 

,250848E-03 

52 

,966451E-01 

,295806E 

01 

603215E 

04 

, 1 0481QE 

08 

,250868E-03 

53 

, 986?QnE-01 

,401 802F 

01 

654949E 

04 

, 1 04837E 

08 

,250892E-03 

54 

,ioo6i3E  no 

,373512E 

01 

707228E 

04 

, 104855E 

08 

,250914E-03 

TEMPFR4TURE 
,SS41AAE  03 
,69U15E  03 
,6l7f.OlE  03 
,S89^S4E  03 
,578983E  03 
.574988E  03 
,574l?4E  03 
.5744S3E  03 
,5ft?884E  03 
,53<599;1E  03 
,5iQ2l8E  03 
*5n0237E  03 
.482872E  03 
.487040E  03 
,4*)2535E  03 
.439272E  03 
.427178E  03 
,41<S177E  03 
.4059S6E  03 
.398133E  03 
.386841E  03 
.378081E  03 
.369828E  03 
,3ft?00QE  03 
,354613E  03 
.34768SE  03 
.341123E  03 
.33489SE  03 
.328979E  03 
.323301E  03 
.317827E  03 
,3l?58lE  03 

.307532E  03 
,3n?786E  03 
.2Q81R3E  03 
.293888E  03 
.2892Q7E  O3 
.28R125E  03 
,28ni7E  03 
.277324E  03 
.273696E  03 
,2701R5F  03 
,28S7R9E  03 
.2634R7E  03 
.260229E  03 
.257092E  03 
,2540ft4E  03 
.251070F  03 
.2481R3E  03 
,245300E  03 
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Table  A-1.  RIP  Output  Summarizing  the  PETROS  3 Input  (Cont'd) 


55 

, 162596E 

00 

.441249E 

01 

J783912E 

04 

, 1 04668E 

08 

.250931E-03 

,242594E 

03 

56 

^in456nE 

00 

,493867E 

01 

,667775E 

04 

, in4676E 

06 

.250943E-03 

.239924E 

03 

57 

.196563E 

00 

,356b67E 

ni 

. 979207E 

04 

, 1 04686E 

06 

.250956E-03 

,237338E 

03 

56 

, 10R546E 

on 

.4695l9t 

ni 

i 1 03625E 

05 

. 104699E 

06 

>250973E-n3 

.234860C 

0 3 

59 

, 1 Ut52‘>E 

no 

.4  34036F. 

ni 

. 10941  7F 

05 

. 1049.1  2E 

06 

,25090nF-03 

.232477F 

03 

60 

,112512E 

00 

.162in9E 

01 

f 1 14404E 

05 

, 104926E 

06 

.251007E-03 

.23ni63E 

03 

5| 

. 1 14495E 

00 

,57n3n6E 

01 

, 1 19429E 

05 

. 104936E 

06 

.25in23E-03 

.227916E 

03 

62 

, 1 16476E 

00 

•164327E 

01 

. 123765E 

05 

.104952E 

08 

. 251 040E-03 

.225715E 

03 

63 

,11646nE 

00 

.715614E 

01 

^ 127665E 

05 

.104965E 

08 

.251057E-03 

.223573E 

03 

64 

,12n442E 

on 

-.153166E 

01 

f 131690E 

05 

, 104976E 

08 

,25l073E-03 

.221453E 

03 

65 

.122425E 

on 

,694621E 

01 

, 134936E 

05 

, 104992E 

06 

.25109nE-03 

.219399E 

03 

66 

• 12«407E 

00 

•266S36E 

01 

, 138159E 

05 

.I05005E 

06 

.251 106E-03 

.217413E 

03 

67 

, 12636QE 

00 

.216212E 

01 

Il4ll76E 

05 

. 105016E 

06 

,251 123E-03 

.215465E 

0 3 

66 

.126371E 

on 

•535449E 

01 

. 145416E 

05 

, 105029E 

nn 

.251 137E-03 

.213551E 

0 3 

69 

, 130353E 

no 

, 15n429E 

nn 

,148743E 

05 

. 105n4OE 

06 

.251 15IE-03 

.21 1702E 

03 

70 

• 132334E 

00 

• 1526S6E 

«i 

-152942E 

05 

. 10505IE 

06 

.251 164E-03 

«2n9864E 

03 

71 

,13a316E 

00 

-.134739E 

01 

, l567bnE 

05 

. 105061E 

06 

.251 177E-03 

.2n8097E 

03 

72 

. 136297E 

00 

«Q24900E 

nn 

’ 1609O4E 

05 

, 105070E 

06 

.251 169E-03 

.2n6367E 

03 

73 

.138279E 

00 

-•451444E 

01 

;i64679E 

05 

. l0506nE 

08 

•251201E-03 

.204649E 

03 

74 

«1 40260E 

00 

-,134501E 

01 

! 168363E 

05 

, 105n89E 

08 

•251213E-03 

.202966E 

03 

75 

• H2241  E 

on 

-.596995E 

01 

. 172018E 

05 

. 105099E 

06 

.251225E-03 

.201347E 

03 

76 

, 1 44222E 

on 

-,422209E 

ni 

; 176166E 

05 

.105107E 

08 

.251236E-03 

. 199720E 

03 

77 

« 146203E 

on 

-,6n2355E 

ni 

, 183667E 

05 

.1051 lOE 

08 

.251241E-03 

« 1 Q61 09E 

03 

76 

, 1 46164E 

on 

-.467372E 

01 

^ 194157E 

05 

, 105106E 

06 

.25124|E-03 

, 196495E 

03 

79 

.150165E 

on 

-.117142E 

01 

,210965E 

05 

.105097E 

08 

.25123nE-03 

.194646E 

03 

60 

, 152146E 

00 

-.166747E 

01 

.228442E 

05 

* l05n85E 

08 

.251216E-03 

.193253E 

03 

61 

.154127E 

00 

-,596696E 

01 

j247503E 

05 

, 105070E 

06 

•25l204E-n3 

.191624E 

03 

62 

, 156106E 

00 

•.267724E 

Ol 

.258163E 

05 

, 1 05066E 

nn 

«25120lE*n3 

•190243E 

03 

63 

.156090E 

on 

••171632E 

01 

.2624n9E 

05 

.105071E 

08 

•251207E*03 

.169013E 

03 

64 

• 160071E 

on 

•,  147660E 

Ol 

.240766E 

05 

•105111E 

06 

.251251E-03 

, 188316E 

03 

65 

, 162052E 

on 

-.305337E 

01 

i?07359E 

05 

.105I63E 

06 

.251309E-03 

, 186219E 

03 

66 

, 164032E 

00 

,453725E 

no 

.1491  12E 

05 

,1Q5105E 

06 

•251326E-03 

.195727E 

03 

• 

THICKNESS  " 

Il61 162E 

on 

spall  fpag 

-HENT  « : 

velocity 

OF 

HASS  center 

■ .267 

167E  03  kinetic 

E ■ .359703E 

04 

J 

X IINCHES) 

TRAtJSVERSF 
STRESS  (PSD 

tangential 

STRESS  IPSD 

bulk 

HOOJLUS  (PSD 

OEnSITT 

TEHPFRATURE 

68 

.165515E  00 

.930598E  01 

;744571E  04 

.105170E  06 

•251295E-03 

.210060E  03 

thickness 

■ .'502340E-03 

SPALL  FRAqHENT  • 

3 

velocity 

OF 

HASS  CENTER 

■ i 

7480n7E  03  1 

KINETIC  C > .366031E 

06 

transverse 

tangential 

SULK 

density 

J 

X unchesi 

STRESS  (PSD 

STRESS  IPSD 

HOOJLUS  IPSD 

TEHPFRATURE 

90 

,l7im6E 

00 

-.632761E 

on 

:323876E 

05 

.105034E 

06 

.251 174E-03 

.182547E 

03 

91 

,173086E 

00 

-.676475E 

Ol 

.248076E 

05 

.105169E 

08 

,251323E-03 

.180676E 

03 

92 

. 175066E 

on 

.257191E 

no 

‘215518E 

05 

.105237E 

06 

.251400E-03 

.178719E 

03 

93 

. 17704RE 

00 

.395966E 

Ol 

.16218IE 

05 

.105331E 

06 

.251504E-03 

.177430E 

03 

94 

,179027E 

on 

-.603212E 

Ol 

‘ 128974E 

05 

. 105391E 

08 

.251572E-03 

.176400E 

03 

95 

.16|nn5E 

nn 

./46244E 

01 

,124377F 

05 

.105411E 

on 

.25159.1E-03 

.175069E 

03 

96 

.1629B3t 

on 

,26928nE 

01 

.130922E 

05 

.105412E 

06 

.251597E-03 

,173811E 

03 

97 

• 1 64961E 

00 

. .4  1799  5F. 

Ol 

I 1 1605SE 

05 

. 105446E 

06 

.251635E-03 

. 1 72566E 

03 

98 

.165939E 

00 

.6024  1 7E 

01 

,665734F 

04 

. 1 05502E 

06 

.251697E-03 

.17M59E 

03 

99 

.1869I6E 

no 

. 157616E 

Ol 

,526292F 

04 

• 105564c 

06 

.25I766E-03 

. 170363E 

03 

1 no 

, l90n93E 

on 

232662F 

01 

,2452?5E 

n4 

« I05M5E 

08 

.25I823E-03 

• 169S12E 

03 

mi 

. l9?87nE 

on 

-.2421  ME 

ni 

-.512S76E 

04 

. 1 n57J5E 

nn 

.251 954E-03 

• 169162E 

03 

m2 

, 1 94645E 

00 

-.2601 99F 

Ol 

-;597528E 

04 

. 105759E 

08 

.251983E-03 

.167659E 

03 

m3 

.196820E 

on 

-.505322E 

nn 

l653nOE 

05 

« 1 05924E 

06 

•252164E-03 

. 167546E 

03 

m4 

•19R794E 

00 

-.306957E 

ni 

-i 193012E 

05 

.105974C 

06 

.252219E-03 

.166773E 

03 

thickness 

■ ^276675E 

-01 
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SUBROUTINE  PDATA  (LINK) 

1 

DIMENSION  TPL0T(50) 

PDAT 

2 

COMMON  (USE  MAIN) 

PDAT 

3 

COMMON  /PDAMP/  NT  I ME ( 50 ) , K I NE T1 ( 20 ) , K I NET2 ( 20 ) , SUMK I N ( 20 ) t 

PDAT 

4 

1CINS(20)»JSS(A0),DST0P(20) »NP , TDAMP , DT ACT » KF , VC EN ( 20 ) »NC0M(20) 

PDAT 

5 

COMMON  /EOSMAT/  GTEMP(IOOO) 

PDAT 

6 

REAL  KINET1.KINET2 

PDAT 

7 

DATA  LTP/0/ 

PDAT 

8 

PDAT 

9 

GOTO  ( I0,20t 50t60t500) .LINK 

PDAT 

10 

10 

READ(5,100)  NP,(NTIME( I),I«1,NP) 

PDAT 

11 

READ  (5,100)  LMAXP 

PDAT 

12 

READ  (5,160)  (TPLOT( I),I=1,LMAXP) 

PDAT 

13 

NT=>1 

PDAT 

14 

KRITE(6,110)  (NTIME( I ) ,I=1,NP) 

PDAT 

15 

WRITE (6, 150)  (TPLOT( I) ,1=1, LMAXP) 

PDAT 

16 

READ  (5,120)  TDAMP, DFACT 

PDAT 

17 

WRITE  (6,130)  TDAMP, DFACT 

PDAT 

18 

NP»1 

PDAT 

19 

KMAX=0 

PDAT 

20 

DO  11  1=1,20 

PDAT 

21 

DST0P(I)=0.0 

PDAT 

22 

KINETK  I )=0.0 

PDAT 

23 

KINET2( I )=0.0 

PDAT 

24 

11 

SUMKIN( I )=0.0 

PDAT- 

25 

NPL0T«3 

PDAT 

26 

PDAT 

27 

IF(NCYCLE  .LE.  DGOTO  45 

PDAT 

28 

12 

CALL  SKIPFILE  (NPLOT.l) 

PDAT 

29 

READ  (NPLOT)  IFLAG 

PDAT 

30 

IF(IFLAG  .EO.  99999)G0T0  14 

PDAT 

31 

READ  (NPLOT)  IFLAG 

PDAT 

32 

IF( IFLAG  .NE.  NCYCLEIGOTO  12 

PDAT 

33 

lA 

CALL  BACKFILE  (NPLOT.l) 

PDAT 

34 

GOTO  60 

PDAT 

35 

PDAT 

36 

20 

CIN1=0.0 

PDAT 

37 

CIN2=0.0 

PDAT 

38 

CIN3=0.0 

PDAT 

39 

CIN4=0.0 

PDAT 

40 

«««**«  COMPUTE  KINETIC  ENERGY  (EACH  REGION)  ♦♦♦♦♦♦♦ 

♦PDAT 

41 

DO  25  1=1,20 

PDAT 

42 

NCOM( I )=o 

PDAT 

43 

25 

CONTINUE 

PDAT 

44 

PDAT 

45 

KF»1 

PDAT 

46 

JK«1 

PDAT 

47 

IFIRST=0 

PDAT 

48 

CINS( 1)=0.0 

PDAT 

49 

VCEN( 1)=0.0 

PDAT 

50 

SMASS=0.0 

PDAT 

51 

KEY=1  , 

PDAT 

52 

PDAT 

53 

DO  35  J=1,MAX2 

PDAT 

54 

IFdVFLAG(J)  .EO.  0)G0T0  35 

PDAT 

55 

HOT  LIQUID  REGION  ( VAPER-L I GUI D ) 

PDAT 

56 

IFdVFLAG(J)  .EG.  33)  C I N 1=C  I N1  ♦.  125*DELT  AM  ( J ) ♦ ( U ( J ) 4U(  J + 1 ) ) **2 

PDAT 

57 

LIQUID  REGION 

PDAT 

58 

IF(IVFLAG(J)  .EG.  32)  C I N2=C I N2+. 1 25*DE LT AM ( J ) ♦ ( U ( J ) +U ( J+ 1 ) ) **2 

PDAT 

59 

MELT  REGION  (LIQUID  SOLID) 

PDAT 

60 

50 


o o 


IF(IVFLAG(J)  •EG.  31)  CIN3=CIN3^, 

125*DELTAM(  J)*  (UU)  ♦UI  Jtl)  ) 

PDAT 

61 

IFdVFLAGU)  .NE.  30)G0T0  35 

PDAT 

62 

IF(  IFIRST  .EG.  0)  JSS(JK)»J 

PDAT 

63 

IF  IP.ST=l 

PDAT 

64 

COMPUTE  KINETIC  ENERGY  OF  EACH 

SOLID  PIECE  (FRAGMENT) 

PDAT 

65 

VELOCITY  OF  MASS  CENTER  OF  EACH  FRAGMENT 

PDAT 

66 

.CINStKF)=CINS(KF)+.125*DELTAK(J)* 

{U(J)+U(J+1) )**2 

PDAT 

67 

SMASS  = SMASS*DELTAMU) 

PDAT 

68 

VCEN(KF)=VCEN(KFH-.5*0ELTAMU)*(U 

( J)»U( J*1 ) ) 

PDAT 

69 

ANY  ZONE  COMBINE 

PDAT 

70 

IFIIBFLAG(J)  .NE.  AIGOTO  28 

PDAT 

71 

NC0M(KF<-1)«1 

PDAT 

72 

JSSUK*!  )”J 

PDAT 

73 

JK=JK*2 

PDAT 

74 

IFIRST=0 

PDAT 

75 

GOTO  35 

PDAT 

76 

PDAT 

77 

?8 

IFIIBFLAGU)  .EQ.  2 )KEY*2 

PDAT 

78 

GOTO  (35,30)tKEY 

PDAT 

79 

30 

VCENIKF ) =VCEN(KF ) /SMASS 

PDAT 

80 

IF(NC0M(KF*1)  .EQ.  OIGOTO  31 

PDAT 

81 

KF=KF*1 

PDAT 

82 

VCEN(KF)=0.0 

PDAT 

83 

31 

JSSIJK^l  )=sJ 

PDAT 

84 

JK*.IK*2 

PDAT 

85 

KF=KF«-l 

PDAT 

86 

KEY  = 1 

PDAT 

87 

IFIKF  .GT.  20)  GOTO  70 

POAT 

88 

CINS(KF)*0.0 

PDAT 

89 

VCEN(KF)=0.0 

POAT 

90 

SHASS=0.0 

POAT 

91 

IFIRST=0 

POAT 

92 

35 

CONTINUE 

POAT 

93 

VCEM(KF)=VCEN(KF) /SMASS 

PDAT 

94 

JSS(JK+1)=MAXZ 

PDAT 

95 

SUM  KINETIC  ENERGY  FOR  KF  FRAGMENTS 

PDAT 

96 

VCENTR=VCEN( 1) 

PDAT 

97 

CINA=0. 

PDAT 

98 

DO  38  J’l.KF 

PDAT 

99 

CINA»CIN4FCINS(J) 

PDAT 

100 

38 

CONTINUE 

PDAT 

lOl 

PDAT 

102 

CINLS=CINA 

PDAT 

103 

CINL=CINAVCIN3 

PDAT 

104 

CINVL  = CINA<-CIN3*CIN2 

PDAT 

105 

CINV=CINA<-CIN3<-CIN2*CIN1 

PDAT 

106 

PDAT 

107 

IFLAG»1 

PDAT 

108 

hRITE(NPLOT)  IFLAG 

PDAT 

109 

WR  ITE(NPLOT)  NC YCLE , TI ME . VCENTR ,KF , ( C I NS ( I ) .I^l.KF) 

PDAT 

110 

POAT 

111 

KHAX=MAXO(KMAXtKF) 

PDAT 

112 

KF'KMAX 

PDAT 

113 

CALL  DAMP 

PDAT 

114 

PDAT 

115 

4A 

IFINCYCLE  .NE.  NT  I ME (NP) ) GOTO  A5 

PDAT 

116 

NP=NP*1 

PDAT 

117 

IFLAG»2 

PDAT 

118 

WRITE(NPLOT)  IFLAG 

PDAT 

119 

WRITEINPLOT)  T I ME . MAXZ , I X U ) , S U ) 

vPU)  «SD(J)  tGTEMPU)  «1VFLAG(J)« 

PDAT 

120 

51 


1J= 1,MAXZ  ) 
A5  RETURN 


C 

50  IFLAG=99999 

WRITE(NPLOT)  IFLAG 
GOTO  ^5 
C 

60  END  FILE  NPLOT 
GOTO 

70  WRITE(6,1A0)  KF,NCYCLE 
IFLAG=99999 
WRITE(NPOLT)  IFLAG 
CALL  FORTRAN 
C 

500  IFILTP  .EQ.  0)G0T0  510 
LTP  = 0 
IFLAG=2 

TIME=TIME-DTNP1 
WR  ITE(NPLOT)  IFLAG 

WR ITE(NPLOT)  T I ME »MAXZ , ( X (J ) , S ( J ) ,P ( J ) , SD ( J ) ,GT EMP ( J ) , I VFL AG ( J ) , 
IJ=1,MAXZ) 

TIME=TIME+DTNPl 
GOTO  A5 
C 

510  IF(NT  .GT.  LMAXPIGOTO  45 

IFITIME  .LT.  TPLOT(NT) )GOTO  45 
DTNPl=TIME-TPLOT(NT) 

TIME=TIME-DTNPl 
LTP=l 
NT=NT+1 
GOTO  45 
C 

100  F0RMAT(16I5) 

110  FORMAT! 5X^33HPL0TS  AT  THE  FOLLOWING  TIME  STEPS/ ( 5X» ( 161 5) ) ) 

120  FORMAT  (2E12.6) 

130  FORMAT  ( 25H  START  DAMPING  AFTER  TI ME t E 1 6 , 5 X , 7HDF ACT  =,E12.6) 
140  FORMAT!*  ENLARGE  ARRAY  CINS  K=  »,I5,»  TIKE  STEP  = *.15) 

150  FORMAT! 5X,23HPL0TS  AT  FOLLOWING  T I ME/5X , ! 5E 1 5. 7 ) ) 

160  FORMAT! 5E15.7) 

END 


PDAT 

121 

PDAT 

122 

PDAT 

123 

PDAT 

124 

PDAT 

125 

PDAT 

126 

PDAT 

127 

PDAT 

128 

PDAT 

129 

POAT 

130 

PDAT 

131 

PDAT 

132 

PDAT 

133 

PDAT 

134 

PDAT 

135 

PDAT 

136 

PDAT 

137 

PDAT 

138 

PDAT 

139 

PDAT 

140 

PDAT 

141 

PDAT 

142 

PDAT 

143 

PDAT 

144 

PDAT 

145 

PDAT 

146 

POAT 

147 

PDAT 

148 

PDAT 

149 

PDAT 

150 

PDAT 

151 

PDAT 

152 

PDAT 

153 

PDAT 

154 

PDAT 

155 

PDAT 

156 

PDAT 

157 

PDAT 

158 

PDAT 

159 

PDAT 

160 
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SUBROUTINE  DAMP 

1 

COMMON  (USE  MAIN) 

DAMP 

2 

COMMON  /PDAMP/  NT  I ME ( 50) , KINETK 20 ) ,KINET2 ( 20 ) , SUHKIN ( 20 ) , 

DAMP 

3 

1CINS( 20) , JSS(AO) ,DST0P(20) ,NP,TDAMPfDFACT,KFfVCEN(20) ,NC0M(20) 

DAMP 

A 

REAL  KINET1,KINET2 

DAMP 

5 

DATA  DCPSI/.C0001A5/ 

DAMP 

6 

DAMP 

7 

DO  5 I=lfKF 

DAMP 

8 

XF(NCOM(I)  .EQ.  1)KINET2(I)=0.0 

DAMP 

9 

KINETl(n=NINET2(I) 

DAMP 

10 

XF(NC0M(I1  .EQ.  DCINSd  )to.O 

DAMP 

11 

KXNET2(I)=CINS(I) 

DAMP 

12 

5 

CONTINUE 

DAMP 

13 

CHECK  FOR  START  OF  DAMPING 

DAMP 

14 

IF1TIME  .LT.  TDAMP)RETURN 

DAMP 

15 

DAMP 

16 

IF(KF  .LT#  2)G0T0  12 

DAMP 

17 

DO  10  I=2»KF 

DAMP 

18 

IF(NCOM(I)  .EQ.  1)VCEN(I)=0,0 

DAMP 

19 

IF(VCEN( I)#EQ#  0*0)G0T0  10 

DAMP 

20 

IF(VCENd-l)  .GT.  VCEN(D)  RETURN 

DAMP 

21 

10 

CONTINUE 

DAMP 

22 

12 

KFINAL=1 

DAMP 

23 

DAMP 

24 

JK  = 1 

DAMP 

25 

DO  30  I=lfKF 

DAMP 

26 

IF(NCOM(  I)  .EQ.  DGOTO  27 

DAMP 

27 

IF<DSTOP(I).NE.  0.0)G0T0  53 

DAMP 

28 

XR(KINET2(T)  .GT.  KINETl ( I) )G0T0  27 

DAMP 

29 

SUNK  I N( I )sSUMKIN( I )+CINS(n 

DAMP 

30 

JS=JSS( JK) 

DAMP 

31 

JSN=JSS( JK+1) 

DAMP 

32 

IF(JS  #EQ*  JSN)GOTO  53 

DAMP 

33 

TRANS=0.0 

DAMP 

34 

DO  52  J=JS,JSN 

DAMP 

35 

TRANS=TRANS+ABS(S( J) ) 

DAMP 

36 

52 

CONTINUE 

DAMP 

37 

TRANS=TRANS/FLOAT( JSN-JS+1) 

DAMP 

38 

TRANS=TRANS*CCPSI 

DAMP 

39 

IF(TRANS  .LE.  DFACT)GOTO  53 

DAMP 

40 

COMPUTE  MOMENTUM  OF  FRAGMENT 

DAMP 

41 

JS=JSS( JK) 

DAMP 

42 

JSN=JSS( JK+1) 

DAMP 

43 

IF(NC0M(I+1)  .NE.  0) JSN=JSS(JK+3) 

DAMP 

44 

DO  25  J=JS»JSN 

DAMP 

45 

UrJ)=VCEN(I) 

DAMP 

46 

25 

CONTINUE 

DAMP 

47 

U(USN+1)=VCEN(I) 

DAMP 

48 

27 

JK=JK+2 

DAMP 

49 

30 

CONTINUE 

DAMP 

50 

RETURN 

DAMP 

51 

DAMP 

52 

53 

DSTOP( I ) =1.0 

DAMP 

53 

IFCKFINAL  .‘GE.  KF)GOTO  55 

DAMP 

54 

KFINAL=KFINAL+1 

DAMP 

55 

GOTO  27 

DAMP 

56 

DAMP 

57 

55 

HRITE(6»1C0)  NCYCLE 

DAMP 

58 

MAXCYL=NCYCLE 

DAMP 

59 

TMAX=TIME 

DAMP 

60 

53 


NTIME(NP)=NCYCLE 
CALL  EDIT 

n 

C 

DAMP 

63 

100  F0RMAT(1H1^10X,29HRUN  SELF-TERMINATED  TIME  STEP, 15) 

DAMP 

6A 

END 

DAMP 

65 

54 


SUERCUTINE  PETiNP 

PETROS-3  INPUT 

COMMON  (USE  MAIN) 

COMMON  /PCAMP/  NTI ME ( 50 ) t KI NE T 1 (20 ) t KI NET2 (20 ) » S U PK I N (20 ) , 
1CINS(20),  JSS(AO)  ,OST0P(20)t,NP,T0AMP,OFACT,KF,VCEN(20  ) ,NC0M(2D) 
COMMON  /ECSMAT/  GTEMPdCOC) 

DATA  0CPSI/.CCCClA5/CMT0IN/.3937/GCT0LB/*93502E-04/ 

WRITE  (6»I0C) 

JK»1 

CO  30  I»1,KF 
J3»JSS( JK ) 

JSN=JSS( JK+I ) 

VCENT»VCEN( I )*CMTCIN 
WRITE(6,110)  I,VCENT,CINS (I ) 

WRITE(6»12C) 

IFIRST*! 

CO  20  J»JSfJSN 

TRANS*S( J )*CCPSI 

TANGS*-(P( J )-0.5*S0( J) )*DCPSI 

TEMX*X( J)*CMTOIN 

BMS»(RHO( J )<C( J )**2)*0CPSI 

RHCP»RHO( J )*GCTCLB 

TEMF*32.C+( 1.8*(GTEMP( J)-273.165) ) 

WRITE! 6, 140)  JtTEMX,TRANS*TANGS»8MS,RH0P,TEMF 
iFdFiRST  .EQ.  1)  THICK1»X(  J) 

IFIRST*0 

THICK»( X( J )-THICKl)*CMTOlN 
20. CONTINUE 
JK«JK42 

IF(JS  .EQ.  JSN)  THICK»(X( JS)-X(JS-l) )*CMTOIN 
WRITE(6,15C)  THICK 
30  CONTINUE 
RETURN 

ICO  FORMAT! 1H1,59X,12HPETR0S  INPUT/) 

110  F0RMAT(//55X,16HSPALL  FRAGMENT  »,I3/ 

I35X» 25HVcLCCITY  OF  MASS  CENTER  » »E1 2. 6 ,1 X , I IHK I N ET I C E *E12.6/) 
130  FORMAT! /5CX, lOHTR AN SVERSE  »5 X ,10HTANGENT I AL 1 8 X EULK , 10 X , 7HCENS 
1/3CX,1HJ,AX,10HX  (INCHES)  , AX, 12HSTRESS  ( P S I ) , 3X  , 1 2HSTRESS  (PSD 
2, 12HM0DULUS  ( P S I ) , 19X , IIHTEMPERATURE ) 

1^0  F0RMAT(28X,I3,6(3X,E12.6) ) 

150  F0RMAT(/52X,11HTHICKNESS  »,E12.6) 

ENC 


PETI 
PET  I 
PET  I 
PETI 
PETI 
PETI 
PET  I 
PETI 
PET  I 
PETI 
PETI 
PETI 
PET  I 
PETI 
PETI 
PET  I 
PETI 
PETI 
PETI 
PETI 
PETI 
PETI 
PETI 
PETI 
PETI 
PETI 
PETI 
PET  I 
PETI 
PETI 
PETI 
PET  I 
PETI 
PET  I 
PET  I 
PETI 
PETI 
ITYPETI 
,3XPETI 
PETI 
PETI 
PETI 
PETI 


6 

7 

E 

9 

1C 

11 

12 

12 

1^ 

15 

le 

17 

18 
19 
2C 
21 
22 
23 
2A 

25 

26 
27 
26 
29 
3C 

31 

32 

33 
3A 

35 

36 

37 
36 
39 
AC 
A1 
A2 
A3 
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APPENDIX  B 


BRL-RIP  PLOTTING  PROGRAM 

The  BRL-RIP  plotting  program  is  independent  of  the  RIP  program.  It 
was  written  to  monitor,  by  plotting,  the  output  needed  in  coupling  RIP 
to  PETROS  3.  The  program  makes  use  of  the  BRLESC  FORTRAN  Plotting  Sub- 
2. 6 

routines  vi?hich  can  easily  be  adapted  to  any  other  plotting  system. 

DESCRIPTION  OF  MAIN  PROGRAM  AND  SUBROUTINES 

The  main  program  reads  the  plotting  data  tape  and  controls  the  flow 
of  information.  If  the  variable  IFLAG  equals  one,  data  is  read  and  stored 
in  the  plotting  data  arrays.  If  it  equals  two,  data  is  read,  converted 
to  proper  units  and  tangential  stress  is  computed,  then  the  program  calls 
subroutine  TPLOT.  When  IFLAG  equals  99999,  the  program  calls  subroutine 
GRAPH . 

If  the  number  of  cycles  is  greater  than  MAXA,  subroutine  GRAPH  is 
called.  If  more  cycles  are  needed,  the  following  data  arrays  must  be 
enlarged:  TIM (MAXA) , E(MAXA),  VCEN(MAXA) , and  PKINE(N) , where  subscript 
N is  equal  to  20  times  MAXA. 

Subroutine  TPLOT  plots  the  profiles  through  the  thickness,  for  the 
transverse  stress,  pressure,  tangential  stress  and  temperature.  See 
Figures  6,  7,  8,  and  9 for  examples  of  the  plotted  output. 

Subroutine  GRAPH  produces  two  graphs.  The  first  is  the  kinetic 
energy  versus  time,  which  is  useful  for  monitoring  the  solution,  for 
the  purpose  of  determining  when  to  start  damping  and  when  to  terminate 
the  run.  The  second  is  the  velocity  of  the  mass  center  versus  time. 

See  Figures  10  and  11  for  examples  of  the  plotted  output. 

Subroutine  PLOTD  will  plot  data  (drawline)  for  every  zone  in  the 
mesh,  but  will  skip  any  spall  zone. 


W.  Coleman  and  J.  V.  Lanaham,  BRLESC  FORTRAN  Plotting  Subroutines," 
Aberdeen  Research  and  Development  Center,  Technical  Report  No.  6, 

June  1970, 
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Description  of  Variables 

Name 

Definition 

B(IOOO) 

Array  used  by  the  BRLESC  FORTRAN  Plotting  Subroutines. 

CINS(IOOO) 

The  array  of  numbers,  corresponding  to  the  kinetic  energy 
of  each  fragment  for  a given  time  step. 

DCPSI 

2 

Conversion  constant  (dynes/cm  to  psi) , 

DX,  DY 

Spacing  between  increments  in  the  X and  Y directions  used 
in  labeling  the  axis. 

GTEMP(IOOO) 

Array  of  numbers  corresponding  to  temperature  profile 
through  the  thickness. 

I 

Index  counter. 

I FLAG 

A flag.  When  it  equals  1,  program  will  read  data  for  the 
timewise  plots.  When  it  equals  2,  the  program  will 
read  data  for  transverse  stress,  pressure,  tangential 
stress,  and  temperature  plots.  When  it  equals  99999, 
the  program  will  plot  the  timewise  plots. 

IVFLAG 

See  Reference  1,  page  507. 

J 

Index  counter. 

KS 

Number  of  spall  fragments  in  the  solid  region. 

LINK 

Index  used  by  a computed  GOTO  statement. 

MAXA 

Size  of  plotting  arrays. 

MAXJ 

Number  of  zones  (See  Reference  1,  page  507). 

MJ 

Index  counter  for  number  of  cycles. 

N 

Index  counter  equal  to  the  maximum  MJ’s. 

NCYCLE 

Counter  on  the  number  of  cycles  (see  Reference  1,  page  506). 

NPLOT 

Magnetic  plotting  tape  input  unit  number. 

PJ(IOOO) 

Array  of  numbers  corresponding  to  pressure  profile  through 
the  thickness. 

E.  H.  Fisher,  G.  A.  Lane  and  E.  A.  Cecil,  "EIF,  A One-Dimensional  Material 
Eesponse  Code  - User 's  Guide, " Systems,  Science  and  Softioare,  Eeport 
3SE-751-I,  September  1972. 
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PKINE(20000) 

The  array  of  numbers,  corresponding  to  kinetic  energy  of 
the  fragments  in  the  solid  region. 

SDJ(IOOO) 

The  array  of  niombers,  corresponding  to  stress  deviator 
through  the  thickness. 

SJ(IOOO) 

The  array  of  numbers,  corresponding  to  stress  through  the 
thickness . 

STEP (2) 

Array  corresponding  to  label  of  current  time. 

SYM.  i = 1,  5 

Arrays  corresponding  to  labels  on  the  axes. 

TEMP 

Temporary  name . 

TIM(IOOO) 

The  array  of  numbers,  corresponding  to  time  history. 

TIME 

Current  value  of  time . 

VCEN(IOOO) 

The  array  of  numbers,  corresponding  to  the  velocity  of 
mass  center  of  the  largest  remaining  solid  region. 

XBAR,  YEAR 

Board  coordinates  of  a point  on  the  plotting  surface. 

XJ(IOOO) 

The  array  of  numbers,  corresponding  to  the  spatial 
coordinates . 

XL 

Length  of  x-axis  in  inches. 

XMAX 

Maximum  x-coordinate . 

XMIN 

Minimum  x-coordinate. 

XPAGE 

Plotting  page  length. 

XS,  YS 

Scales  in  the  X and  .Y  directions. 

XSCALE 

Scale  in  the  X direction. 

x^,  X2 

Plotting  arrays  for  the  X and  Y directions. 
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c 

BRL-RIP  PLOTTING  PROGRAM  {BRLESC  FORTRAN  PLOTTING  SUBROUTINES) 

MAIN 

1 

COMMON  TIH( ICOO) tGTEMPI 1000) t I VF LAG  II 000) t VCEN 11000 ) t 

MAIN 

2 

IX J ( 1000)  ,SJ{  1000)  .PJUOOO)  tSDJ{  1000)  ,MAXZ,HJ,NCYCLEf  TIKE 

MAIN 

3 

COMMON  K S , MA X A , C I N S ( 20 ) ♦ PK I ME ( 20000 ) 

MAIN 

4 

DATA  DCPSI/.CCC01A5/ 

MAIN 

5 

♦ 

MVPRTO(OIKOO) 

MAIN 

6 

♦ 

SUBR(PLTCCB^OIAOO) 

MAIN 

7 

c 

MAIN 

8 

NPL0T=3 

MAIN 

9 

REWIND  NPLOT 

MAIN 

10 

MJ=0 

MAIN 

11 

MAXA=1000 

MAIN 

12 

10 

READ(NPLOT)  IFLAG 

MAIN 

13 

IFdFLAG  .EQ.  999Q9)G0T0  AO 

MAIN 

lA 

G0T0(20»25)»IFLAG 

MAIN 

15 

c 

MAIN 

16 

20 

MJ=MJ+1 

MAIN 

17 

RE  AD (NPLOT)  NC YCLE »T I M ( M J ) , VCEN ( MJ ) , KS , ( C I NS ( I) t 1= 1 , KS ) 

MAIN 

18 

IF(NCYCLE  ,GE.  MAXA)  CALL  GRAPH 

MAIN 

19 

c 

MAIN 

20 

TEMP=0.0 

MAIN 

21 

DO  22  I=l»20 

flAIN 

22 

J*MJ^MAXA+( I-l ) 

MAIN 

23 

TEMP=TEMP+CINS( I ) 

MAIN 

2A 

PKINE( J)=TEMP 

MAIN 

25 

22 

CONTINUE 

MAIN 

26 

GOTO  10 

MAIN 

27 

c 

MAIN 

28 

25 

READ  (NPLOT)  T I ME » MAXZ » ( X J ( I ) , S J ( I ) , P J ( I ) » SD J ( I ) » GTEMP ( I ) . 

MAIN 

29 

1IVFLAG(  I ),I  = 1»MAXZ) 

MAIN 

30 

TIME=TIM(MJ) 

MAIN 

31 

DO  30  I=l»MAXZ 

MAIN 

32 

SJ  ( I )=SJ( I )«CCPSI 

MAIN 

33 

PJ ( I )=PJ( I )«DCPSI 

MAIN 

3A 

SDJ( I )=SDJ( I )*DCPSI 

MAIN 

35 

SDJ( I)=PJ( I)-0.5«SDJ(I) 

MAIN 

36 

30 

CONTINUE 

MAIN 

37 

CALL  TPLCT 

MAIN 

36 

GOTO  10 

MAIN 

39 

c 

MAIN 

40 

AO 

CALL  GRAPH 

MAIN 

41 

CALL  EXIT 

MAIN 

42 

END 

MAIN 

43 

60 


^ rg  m -4* 


SUBROUTINE  TPLOT 

DIMENSION  B( ICOO ) , SYMl ( 3 ) t S YM2 ( 3) , SYM3 ( 2 ) t STEP C 2 ) t SYM4 ( 3 ) t SYM5 (2  ) 
COMMON  (USE  MAIN) 

COMMON  /PLOT/B 


DATA  (SYMKI) 
CSYM2( I) 
ISYM3( I ) 
l$YMA( I ) 
C$YM5( I ) 
XL*7.0 
YL*6.0 
L1NK>1 


I*lf  3)/lOHSPACIAL  COf lOHORDINATE  (fAHCM)>/, 
I«lt3)/lOHTRANSVER$E,lOH  STRESS  (P,AHSI)>/f 
I»l,2)/lOHPRESSURE  (,5HPSI)>/ 

I*lf 3)/lOHTANGENTIALtlOH  STRESS  (PfAHSI)>/t 
I«lf 2)/lOHTEMPERATUR|2HE>/ 


XPAGE»12.0 
CALL  PLTCCB 
N«MAXZ 


(XPAGEf IvB(l) .B( 1000) ) 
GRAPH  ONE 


XBARs 

YBAR» 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 


3.0 

2.0 

FIXSCA 

FIXSCA 

PLTCCS 

PLTCCA 

PLOTD  ( 

LABELA 

PLTCCS 

PLTCCT 

PLTCCT 

PLTCCT 


(XJ( I ) fNfXLf  XStXMINtXMAXfOX) 
(SJ( 1) tN,YLtYStYMIN,YMAX,DY) 
(XBARtYBAR ,XM I N, YM I N , XS t YS ) 
(OX,DY,XMINtXMAXf YMINtYMAXf A) 
XA . %A  \ 


50  CALL  PLTCCT  ( . I t 6HCYCL E> . 0 . 0 , I . 
’ ENCODE  (6,100, STEP)  NCYCLE 
CALL  PLTCCT  ( . I , STEP ( 1 ) , 0. 0 , I . 0 
CALL  PLTCCT  ( . 1 , 6HT I ME=>, 0. 0 , 1. 
ENCOOE( 12, llOtSTEP)  TIME 
CALL  PLTCCT  ( . I , STEP ( I ) , 0.0, 1 . 0 
GOTO( 55,60,65,70 .LINK 

GRAPH 

55  YBARST12.0 

CALL  FIXSCA  (PJ(l),N,YL,YS,YMIN 
CALL  PLTCCS  ( XBAR , YBAR , XM IN , YM I 
CALL  PLTCCA  ( OX , OY , XMI N, XMAX , YM 
CALL  PLOTO  (XJ,PJ) 

CALL  LABELA  ( OX , OY , XMI N, XMAX , YM 
CALL  PLTCCS  ( XBAR , YBAR , 0. 0, C . 0 , 
CALL  PLTCCT  ( . I , SYM I ( I ) , 0. 0 , I , 0 
CALL  PLTCCT  ( . I , S YM3 ( I ) , I . 0, 0 . 0 


0,0.0, -1.0) 

,0.6,~1.0) 
0,0. 0,-1. 2) 

,0.6,— 1.2) 

TWO 


,YMAX,OY) 

N,XS,YS) 

IN,YMAX,A) 

IN,YMAX,1.0,1.0) 

1.0, 1.0) 

,2, 3, -0.6) 

,-1.2,2. 3) 


LINK=2 
GOTO  50 


GRAPH  THREE 


60  YBAR 

= 22.0 

CALL 

FIXSCA 

CALL 

PLTCCS 

CALL 

PLTCCA 

CALL 

PLOTO  ( 

CALL 

LABELA 

CALL 

PLTCCS 

CALL 

PLTCCT 

CALL 

PLTCCT 

CALL 

PLTCCT 

[N,YMAX,OY) 

[N,XS,YS) 

UM,YMAX,A) 

UM,YMAX,l.O,  l.O) 

a. 0,1.0) 

),2.3,-0.6) 

1.0,0. 0,-1. A, 2. 5) 
),-1.2, l.R) 


♦ ♦♦♦ 

1 TPLO 

TPLO 

3 

TPLO 

A 

TPLO 

5 

TPLO 

6 

TPLO 

7 

TPLO 

B 

TPLO 

9 

TPLO 

10 

TPLO 

11 

TPLO 

12 

TPLO 

13 

TPLO 

lA 

TPLO 

15 

TPLO 

16 

17 

TPLO 

18 

TPLO 

19 

TPLO 

20 

TPLO 

21 

TPLO 

22 

TPLO 

23 

TPLO 

2A 

TPLO 

25 

TPLO 

26 

TPLO 

27 

TPLO 

26 

TPLO 

29 

- TPLO 

30 

TPLO 

31 

TPLO 

32 

TPLO 

33 

TPLO 

3A 

TPLO 

35 

TPLO 

36 

TPLO 

37 

—TPLO 

38 

TPLO 

39 

TPLO 

AO 

TPLO 

A1 

TPLO 

A2 

TPLO 

A3 

TPLO 

AA 

TPLO 

A5 

TPLO 

A6 

TPLO 

A7 

TPLO 

AB 

TPLO 

A9 

— TPLO 

50 

TPLO 

51 

TPLO 

52 

TPLO 

53 

TPLO 

5A 

TPLO 

55 

TPLO 

56 

TPLO 

57 

TPLO 

58 

TPLO 

59 

TPLO 

60 

61 


LINK>=3 
GOTO  50 


C ' GRAPH  FOUR 

65  CALL  PLTCCP 
XBAR=3.0 
YBAR=2.0 

CALL  FIXSCA  ( GTEHP ( 1 ) , N, YL t YS , YMI N, YM AX , DY ) 

CALL  PLTCCS  ( XB AR , YB AR , XM IN , YP I M , XS , YS ) 

CALL  PLTCCA  ( DX  , D Y , XMI M, XMAX t YH I N , YMAX 1 4 ) 

CALL  PLOTD  (XJ.GTEMP) 

CALL  LABELA-  ( DX , DY , XHI N, XMAX , YM I N, YMAX , 1 . 0, I . 0 ) 
CALL  PLTCCS  ( XB AR , YB AR , 0. C,  0. 0 , 1 . 0, 1 . 0) 

CALL  PLTCCT  ( . 1 , S YM I ( 1 ) , 0. 0 , 1 , 0 , 2. 3 , -0. 6 ) 

CALL  PLTCCT  ( . 1 , SYH5 ( 1 ) , 1 . 0, 0.0 , - 1 . 0, 2. 5 ) 

LINK=A 
GOTO  50 
70  CALL  PLTCCP 
RETURN 
C 

100  FORMAT! I A.2H  >) 

110  F0RMAT(E12.7,2H  >) 

END 


m 

n 

63 

TPLO 

6A 

TPLO 

65 

TPLO 

66 

TPLO 

67 

TPLO 

68 

TPLO 

69 

TPLO 

70 

TPLO 

71 

TPLO 

72 

TPLO 

73 

TPLO 

7A 

TPLO 

75 

TPLO 

76 

TPLO 

77 

TPLO 

7B 

TPLO 

79 

TPLO 

ao 

TPLO 

81 

TPLO 

82 
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SUBROUTINE  GRAPH 

1 

DIMENSION  B( 1000) tSYMl(2) tSYH2(2J,SYM3(AJ 

GRAP 

2 

COMMON  (USE 

MAIN) 

GRAP 

3 

COMMON  /PLOT/B 

GRAP 

4 

DATA  (SYMK  1 )t I=1|2)/10HTIME  ( NAM0t9HSEC0NDS ) >/ t 

GRAP 

5 

1 (SYM2( n, I=lt2)/10HKINETIC  EN,5HERCY>/ 

GRAP 

6 

DATA  (SYM3( I ), I=1,A)/10HVELCCITY  O.IOHF  HASS  CEN.IOHTER  (CM/SE 

♦3HCGRAP 

7 

!)>/ 

GRAP 

8 

GRAP 

9 

N*HJ 

GRAP 

10 

XL«7.0 

GRAP 

11 

YL»6.0 

GRAP 

12 

XSCALE=1,0E9 

GRAP 

13 

1 A. 

XBAR=3.0 

GRAP 

15 

YBAR=2.0 

GRAP 

16 

CALL  FIXSCA 

( TIM( D.N.XL.XS.XHIN.XHAX.DX) 

GRAP 

17 

CALL  PLTCCS 

( XBAR,YBAR,0.0,0.0tl.0tl.0) 

GRAP 

16 

CALL  PLTCCT 

( .ltSYMI(l),0.0,1.0,2.5t-0.6) 

GRAP 

19 

CALL  PLTCCT 

( .I,SYH2( 1),1.0,0.0,-1.2,2.3) 

GRAP 

20 

II 

o 

o 

Q 

GRAP 

21 

J=  UMAXA*(  I- 

1) 

GRAP 

22 

IF(J  .EO*  1) 

CALL  FIXSCA  (PKINE(J),N,YLtYS,YMIN,YMAX,DY) 

GRAP 

23 

CALL  CONSCA 

(PKINE( J) tN,YL,YS,YHIN,YHAX,DY) 

GRAP 

24 

10  CONTINUE 

GRAP 

25 

GRAP 

26 

CALL  PLTCCS 

(XBAR,YBAR,XHIN,YMIN,XS,YS) 

GRAP 

27 

CALL  PLTCCA 

(DX,DY,XMINtXMAXtYMINtYMAXtA) 

GRAP 

28 

DO  ?0  I*UKS 

GRAP 

29 

Js 1+MAXA»( I- 

1) 

GRAP 

30 

CALL  PLTCCD 

( ltO,TlM( 1) ,PKINE( J) ,N,OtXHIN,XHAX,YMlN,YHAX) 

GRAP 

31 

20  CONTINUE 

GRAP 

32 

CALL  LABELA 

(DX ,DY,XHIN,XHAX,YHIN,YHAX,XSCALE,1 .0) 

GRAP 

33 

Y8AR=12.0 

GRAP 

35 

CALL  PLTCCS 

(XBAR,YBAR,0.0,0.0, I.O.I.O) 

GRAP 

36 

CALL  PLTCCT 

( . I.SYMK  I)  ,0.0,  I.0.2.5.-0.6) 

GRAP 

37 

CALL  PLTCCT 

( . I,SYH3( 1 ),I.0,0.0,-I.2,1.3) 

GRAP 

36 

CALL  FIXSCA 

( VCEN( 1),N,YL,YS,YHIN,YMAX,0Y) 

GRAP 

39 

CALL  PLTCCS 

( XBAR , YBAR,XHIN,YHIN,XS,YS) 

GRAP 

40 

CALL  PLTCCA 

( DX,DY,XHIN, XHAX,YHIN,YMAX,A) 

GRAP 

41 

CALL  PLTCCD 

(1,0,TIH(1),VCEN(1),N) 

GRAP 

42 

CALL  LABELA 

(DX,DY,XHIN,XHAX,YMIN,YMAX,XSCALE,1.0) 

GRAP 

43 

CALL  PLTCCP 

GRAP 

44 

RETURN 

GRAP 

45 

END 

' 

GRAP 

46 
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SUBROUTINE  PLOTO  (Y1»Y2) 

1 

COMMON  (USE  MAIN) 

PLOT 

2 

DIMENSIQN  XI ( 1000) ,X2( 1000) ,Yl( 1) tY2(l) 

PLOT 

3 

PLOT 

4 

K*0 

PLOT 

5 

DO  5 J=l»MAX2 

PLOT 

6 

IFUVFLAG(J)  .EO.  0)G0T0  10 

PLOT 

7 

K*K41 

PLOT 

8 

Xl(K)=Yl(J) 

PLOT 

9 

X2(K)=Y2(J) 

PLOT 

10 

5 

CONTINUE 

PLOT 

11 

10 

CALL  PLTCCD  ( 1 , 0, X 1 ( 1 ) ,X2 ( 1 ) t K ) 

PLOT 

12 

K*0 

PLOT 

13 

IF(J  .GE.  MAX2)  RETURN 

PLOT 

GOTO  5 

PLOT 

15 

END 

PLOT 

16 
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APPENDIX  C 


CODING  CHANGES  TO  PETROS  3 

Two  changes  were  required  in  the  BRL  version  of  PETROS  3 for  it  to 
be  able  to  accept  RIP  data.  The  first  change  was  the  alteration  of  the 
material  properties  subroutine,  MATPRO,  to  permit  the  introduction  of 
variable  material  properties  (E,  v,  p,  a ) for  each  Gaussian  layer  through 
the  thickness.  ^ 

The  second  change  was  made  in  order  to  have  an  initial  stress 
distribution  through  the  thickness  (see  Figure  8)  rather  than  a zero 
stress  state  at  the  start  of  the  calculation.  The  subroutine  ISTRES 
reads  spatial  coordinates  and  tangential  stress  (see  Table  A-1) . After 
all  cards  are  read,  linear  interpolation  through  the  thickness  for  values 
of  tangential  stress  at  each  Gauss  point  is  performed.  Then  the  initial 
strain  energy  is  determined.  This  subroutine  is  called  from  within 
subroutine  STEPl,  following  card  STP11020,  page  349,  Reference  5,  by  the 
statement  CALL  ISTRES;  it  is  only  called  at  time  step  zero. 

The  following  is  a FORTRAN  listing  of  the  two  subroutines  as  they 
are  used  in  the  BRL  version  of  PETROS  3. 
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u u u 


SUBROUTINE  MATPRO 
EVALUATE  MATERIAL  PROPERTIES 
THIS  VERSION  IS  FOR  RIP  MATERIAL  PROPERTIES 

COMMON  (USE  MAIN) 

COMMON  /CONTN2/  EE  1 ( A ) , E E? ( A 1 , HNUl ( A 1 , HNU2 ( A ) , 

lCONSTZ(Al,CONSTI( A),C0NST2(A1  .tXPONZ(A)  .EXPONl ( A 1 , E XPON2 ( A 1 , 

2 RHOl ( A1 ,RH02( A) , SIG1DZ( A,5 ) .SIGUU  ( A.5)  , S I Gl 02 ( A, 5 ) , 

3EPS10Z(A,51 ,EPS1D1(A,5),EPS1D2(A,5),C0EFFZ(A,51,SIGMAZ(A,51, 
ATHIKNS( A,  1A,21  1 ,CENTRZ( A, IA.2I 1 .ALPHAZ ( A)  ,E£Z (6  ,6  1 ,HNUZ (5,61  , 
5RH0Z(6,61 

DIMENSION  SIGIDTI  51,EPS10T(  51,EEET(  61 
C 

IFIOMATPRl  GC  TO  103 
WR  ITEINWR  ITE, 10  1 

10  FORMATI'l  SUBROUTINE  MATPRO  RIP  MATERIAL  PROPERTIES'! 

11  FORMAT! 5£15. 61 

00  101  ILAY=  I.NLAYER 

READ  (MREAO.Ill  CCNSTZ ( I LAYl ,EXP0NZ ( I LAYl , ALPHAZ ( I L AYl 
WRITE  (N«RITE,1I1  C ON  S TZ  ( I L A Y 1 , EXPONZ  ( 1 L A Y 1 , AL  PH  AZ  ( I L AY  l' 
NGAUSL=MGAUSS( ILAYl 
00  101  IGAU  =I,NGAUSL 

READ  (NREAD.lll  RHOZ ( I LA Y , I G AU  1 , EFZ ( 1 L A Y , I G AU  1, 

1HNUZ( ILAY, IGAU  1 , S I Gl DZ ( I L A Y , I GAU  1 
WRITE  (NWRITE.lIl  RHOZ ( ILAY, IGAU  1 , EE Z ( I L A Y , I G AU  1, 

1HNUZ(  ILAY, IGAU  1 , S I G 1 DZ ( I L A Y , I G AU  1 
ICl  CONTINUE 

(3MATPR=.  TRUE. 

C 

103  RHC=RHOZ ( ILAYER, IGAUSSl 
EE=£EZ(  ILAYER, IGAUSSl 
HNU  = HNUZ(  ILAYER, I GAUSS  1 
CONST=CONSTZ ( ILAYER) 

EXPON=EXPONZ ( ILAYER) 

ALPHA=ALPHAZ (ILAYER 1 
C 

NSBL=NSUQL( ILAYER 1 
DO  100  ISUBL=1,NSBL 
COEFFI  ISUDL 1 = 1.0 

SI GHA(  ISU3L )=SIG1DZ( ILAYER, IGAUSSl 
100  CONTINUE 
RETURN 
END 
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MATP 
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MATP 
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MATP 
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MATP 
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MATP 
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MATP 

1 A 

MATP 
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MATP 

16 

MATP 
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MATP 

MATP 
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MATP 
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MATP 
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MATP 
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MATP 
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MATP 
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MATP 
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MATP 
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MATP 
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MATP 

26 

MATP 
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MATP 

30 

MATP 

31 

MATP 

32 

MATP 

33 

MATP 

34 

MATP 

35 

MATP 

36 

MATP 

37 

MATP 

33 

MATP 

39 

MATP 

40 

MATP 

41 

MATP 

42 

MATP 

43 
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SUBROUTINE  ISTRGS 

1 

COMMON  (USE  MAIN) 

ISTRE 

2 

COMMON  /CCNTN3/  TAUll(lAt2lt  6 ) , TAUl 2 ( 1 A , 2 1 , 6 ) , TAU22 ( 1 A t 2 1 1 6) 

ISTRE 

3 

DIMENSION  X(  ICC) , STRES( 100) 

ISTRE 

ISTRE 

A 

5 

IF (QSTRES)  GCTO  3C 

ISTRE 

6 

QSTRES=«TRUE. 

ISTRE 

7 

READ  INITIAL  STRESS  DATA  FROM  CARDS 

ISTRE 

8 

READ  (5f5CO)  N 

ISTRE 

9 

DO  5 1=1, N 

ISTRE 

10 

REAO  (NREA0,510)  X ( I ) , STRE S ( I ) , OUM 

ISTRE 

11 

X(  I ) = X( I )-X(  1) 

ISTRE 

12 

5 

CONTINUE 

ISTRE 

13 

HALF=.5*X(N ) 

ISTRE 

lA 

NN=N“1 

ISTRE 

15 

X(1)=HALF 

ISTRE 

16 

X(N)=-HALF 

ISTRE 

17 

DO  10  1=2, NN 

ISTRE 

la 

X(  I )=HALF-X(  I) 

ISTRE 

19 

10 

CONTINUE 

ISTRE 

20 

ISTRE 

21 

WRITE(6,520) 

ISTRE 

22 

DO  20  1=1, N 

ISTRE 

23 

WRITEI6, 530)  X( I ) ,STRES( I ) 

ISTRE 

2A 

20 

CONTINUE 

ISTRE 

25 

ISTRE 

26 

30 

IF(Z.LT.X(N)  .OR,  Z.GT.X( 1) )GCTO  50 

ISTRE 

27 

CO  AO  1=2, N 

ISTRE 

28 

IF(Z  .GE.  X(  I)  )GOTO  A5 

ISTRE 

29 

AO 

CONTINUE 

ISTRE 

30 

A5 

IM*I-l 

ISTRE 

31 

FSTRE=STRES( IM)+(Z-X( IM) ) ♦ ( STRE S ( I ) -STRE S ( I M ) ) / ( X ( I ) -X ( IM ) ) 

ISTRE 

32 

GOTO  55 

ISTRE 

■a  -a 

50 

FSTRE=O.C 

ISTRE 

3A 

55 

00  151  LA=1,2 

ISTRE 

35 

CO  151  LB=1,2 

ISTRE 

TS (LA,LB  ) = C.O 

ISTRF 

37 

151 

G(LA,LB)=A(LA,LB)-2.«Z«B(LA,LR) 

ISTRE 

38 

GO  = C(  1,1 )*G( 2,2)-G( 1,2)**2 

ISTRE 

39 

SRG=SQRT( GO) 

ISTRE 

ISTRE 

^0 

A1 

NSEL=NSUBL( ILAYER ) 

ISTRE 

A2 

CO  160  ISB=1,NSRL 

ISTRE 

A3 

IZ  = IZ  + 1 

ISTRE 

AA 

TAUIK 11, 12, IZ)=FSTRE/G(1,  1) 

ISTRE 

A5 

TAU12( 11, 12, IZ)=0,0 

ISTRE 

A6 

TAU22(  11,  12,  IZ)=FSTR£/G(2,2) 

ISTRF 

A7 

TAUSBU  1,  1)  = TAU11(  I1,I2,IZ) 

ISTRE 

AB 

TAUSDK 1,2)=TAU12( 11, 12, IZ) 

ISTRE 

A9 

TAU3BU2, 2)  = TAU2?(  11,  12,  I Z) 

ISTRE 

50 

TAUSBL(2,1)=TAUSBL(1,2) 

ISTRE 

ISTRE 

f 1 
52 

DO  155  L A = l, 2 

ISTRE 

55 

CO  155  LB=1,2 

ISTRE 

5A 

TM(LA,LB)=TALSBL(LA,  l)*r.(  1 , LB  ) ^TAUSBL  (LA  , 2 ) ( 2 , LB  ) 

ISTRE 

c ^ 

153 

TS(LA,LR)=TS(LA,LB)  + COE FF ( I SB ) ♦ TM ( L A , LP ) 

ISTRE 

5 6 

160 

CO NT INUE 

ISTRE 

*^7 

strain  ENERGY 

ISTRE 

58 

NGAUSL=NGAUSS(  ILAYER  ) 

ISTRE 

59 

67 


PAR=THICKN^. 5*WEIGHT( IGAUSS,NGAUSU) 

CB  = PAR*CA/EE 

IFdl.EO.ISTRl  .OR.  I1.E0.ISTR3)SRG  = 0.5*SRG 
IFn2.ECdSTR2  .OR.  I2.PU.  I STRA  ) SRC  = 0.  5*SRG 

STREN  = STREN4(  (TSn,l)  + TS(2f2)  )^*2~2.*(  1.0>HNU)^(TS(  1 f 1)  «TS(2t2) 
I - TS( 1, 2)^TS(2, 1) ) )^SRG*CB 

RETURN 


500  FORMAK  16) 

510  FORMAK  3E15. 6) 

520  FORMAK IHl, 52X, 28HSTRESS  THROUGH 
153X,9HTHICKNESS,  ICX, 6HSTRESS ) 

530  F0RMAT(51X,E13.6,AX,E13.6) 

ENC 


THE  THICKNESS/70X,10HTANGENTAIL/ 


ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 

ISTRE 
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APPENDIX  D 


TEMPERATURE  VARIATION  DURING  STRUCTURAL  RESPONSE 

As  noted  previously,  the  PETROS  3 calculations  were  made  on  the 
premise  that  the  initial  temperature  distribution  through  the  thickness 
of  the  residual  plate  was  frozen  for  the  duration  of  the  structural  re- 
sponse prediction  (5,520  microseconds).  Although  this  duration  is  quite 
brief,  it  was  felt  desirable  to  determine  the  extent  of  temperature  re- 
distribution which  could  occur  during  this  period.  Consequently,  a heat 
transfer  problem  for  the  residual  aluminum  plate  was  formulated.  By 
assuming  that  the  edges  of  the  plate  were  insulated,  the  problem  is  re- 
duced to  one  spatial  dimension:  the  through-thickness  direction.  On 

the  front  and  back  faces  of  the  residual  plate  the  boundary  conditions 
were  taken  to  be  the  radiation-type.  The  solution  to  this  intial-value 

1 7 

problem  was  obtained  by  using  the  CINDA  code  . The  results,  consisting 
of  a comparison  of  temperature  distributions  at  the  initial  and  final 
times  of  the  PETROS  3 calculation,  are  shown  in  Figure  D-1.  The  abscissa 
in  this  plot  represents  the  distance  from  the  front  face  of  the  residual 
plate.  It  may  be  seen  that  there  is  an  appreciable  change  in  temperature 
during  the  plate  response,  especially  in  the  neighborhood  of  the  front 
face . 


It  would  obviously  be  desirable  to  take  this  temperature  variation 
into  account  in  future  structural  response  calculations  of  this  type. 

This  may  readily  be  accomplished  if  one  is  willing  to  make  the  plausible 
assumption  that  thermomechanical  coupling  can  be  neglected.  The  CINDA 
code  could  be  employed  to  obtain  the  transient  temperature  distributions 
(at  frequent  time  intervals)  and  storage  of  these  data  on  a magnetic 
tape  could  be  read  by  the  PETROS  3 code  as  the  structural  response  solu- 
tion is  marched  out  in  time. 

Recalling  that  the  PETROS  3 calculations  were  made  using  both  cold 
and  heated  material  properties  in  an  effort  to  bound  the  true  solution, 
the  latter  case  corresponding  to  instant  degradation  of  material  strength 
with  temperature,  raises  another  question  as  to  the  proper  modeling  of 
the  physical  event  when  heat  transfer  is  considered.  That  is,  should 
the  instantaneous  recovery  of  material  strength  be  allowed  in  those  por- 
tions of  the  structure  which  are  cooking?  Full  recovery  of  strength  does 
not  occur  in  metals  which  are  subjected  to  elevated  temperatures  for  long 
soak  times  and  then  slowly  permitted  to  cool.  While  the  response  under 
consideration  certainly  does  not  allow  time  for  annealing,  it  would  be 
more  consistent  with  the  objective  of  predicting  a bound  on  the  response 
in  the  heated  material  case  if  one  allowed  no  recovery  of  material  strength 
following  the  maximum  temperature  excursion. 
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J,  D.  Gask^  "Chrysler  Improved  Numerical  Differencing  Analyser  for 
Third  Generation  Computers^"  Chrysler  Corporation  Space  Division^ 
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